\LargeHardy-Schatten Norms of Systems, Output Energy Cumulants and Linear Quadro-Quartic Gaussian Control
Igor G. Vladimirov
\qquad Ian R. Petersen \thanks{This work is supported by the Australian Research Council. The authors are with the School of Engineering and Information Technology
Canberra ACT 2600
This paper is concerned with linear stochastic control systems in state space. The integral of the squared norm of the system output over a bounded time interval is interpreted as energy. The cumulants of the output energy in the infinite-horizon limit are related to Schatten norms of the system in the Hardy space of transfer functions and the risk-sensitive performance index. We employ a novel performance criterion which seeks to minimize a combination of the average value and the variance of the output energy of the system per unit time. The resulting linear quadro-quartic Gaussian control problem involves the $\cH_2$ and $\cH_4$-norms of the closed-loop system. We obtain equations for the optimal controller and outline a homotopy method which reduces the solution of the problem to the numerical integration of a differential equation initialized by the standard linear quadratic Gaussian controller.
quadro_quartic_MTNS_arxiv.tex
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%2345678901234567890123456789012345678901234567890123456789012345678901234567890
%        1         2         3         4         5         6         7         8

\documentclass[letterpaper,12pt,conference]{ieeeconf}  % Comment this line out
                                                          % if you need a4paper
%\documentclass[a4paper, 10pt, conference]{ieeeconf}      % Use this line for a4
                                                          % paper

\IEEEoverridecommandlockouts                              % This command is only
                                                          % needed if you want to
                                                          % use the \thanks command
\overrideIEEEmargins
% See the \addtolength command later in the file to balance the column lengths
% on the last page of the document

%Page Settings
%  	Letter Paper 	A4 Paper
%Paper size	(8.5x11.0)in 	(21.59x27.94)cm 	(8.27x11.69)in 	(21x29.7)cm
%Top margin (1st page) 	1.0in	2.54cm 	1.0in	2.54cm
%Top margin (rest) 	.75in 	1.9cm 	0.75in 	1.9cm
%Left margin 	0.75in	1.9cm 	0.75in	1.9cm
%Right margin 	0.75in	1.9cm 	0.52in	1.32cm
%Bottom margin 	0.7in	1.8cm 	1.44in 	3.67cm
%Text width 	7.0in 	17.78cm 	7.0in 	17.78cm
%Text height 	9.5in	24.13cm 	9.5in	24.13cm
%Column width 	3.4in 	8.63cm 	3.4in 	8.63cm
%Column separation 	0.2in 	0.5cm 	0.2in 	0.5cm
%	  	
%Font Settings
%Title	16pt
%Author list	11pt
%Section titles	10pt
%Normal text	10pt
%Table/Figure captions	10pt
%Sub-scripts	8pt
%Sub-sub-scripts	7pt
%Footnotes	8pt



% The following packages can be found on http:\\www.ctan.org
%\usepackage{graphics} % for pdf, bitmapped graphics files
%\usepackage{epsfig} % for postscript graphics files
%\usepackage{mathptmx} % assumes new font selection scheme installed
\usepackage{times} % assumes new font selection scheme installed
%\usepackage{amsmath} % assumes amsmath package installed
%\usepackage{amssymb}  % assumes amsmath package installed


\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{graphicx}
\usepackage{color}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\setlength{\textwidth}{175truemm}
%\setlength{\oddsidemargin}{-6.5truemm}
%\setlength{\evensidemargin}{1.6mm}
%\setlength{\textheight}{260truemm}
%
%\setlength{\topmargin}{-25truemm}


%\setlength{\textwidth}{170truemm}
%\setlength{\oddsidemargin}{-5truemm}
%\setlength{\textheight}{250truemm} \setlength{\topmargin}{-20truemm}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newtheorem{definition}{Definition}
\newtheorem{lemma}{Lemma}
\newtheorem{theorem}{Theorem}
\newtheorem{proposition}{Proposition}
\newtheorem{corollary}{Corollary}
\newtheorem{heuristic}{Heuristic}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\<{\leqslant}           % nice less than or equal to sign
\def\>{\geqslant}           % nice larger than or equal to sign
\def\div{{\rm div}}         % divergence
%\def\endproof{{$\Box$}}
\def\endremark{{$\blacktriangle$\bigskip}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\d{\partial}
\def\wh{\widehat}
\def\wt{\widetilde}


\def\Re{{\rm Re\,}}   % real part
\def\Im{{\rm Im\,}}   % imaginary part

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\cH{{\cal H}}   % Hardy space

\def\mZ{{\mathbb Z}}    % set of integers
\def\mN{{\mathbb N}}    % set of positive integers
\def\mR{{\mathbb R}}    % real line
\def\mC{{\mathbb C}}    % complex plane

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\llaw{{[\![}}       % probability law
\def\rlaw{{]\!]}}       % probability law
\def\Tr{{\rm Tr}}       % matrix trace
\def\rT{{\rm T}}        % matrix transpose
\def\rF{{\rm F}}        % matrix transpose

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\def\bP{{\mathbf P}}    % probability
\def\bE{{\mathbf E}}    % expectation

\def\bK{{\mathbf K}}    % cumulant

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\def\bra{{\langle}}
\def\ket{{\rangle}}

\def\Bra{\left\langle}
\def\Ket{\right\rangle}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\def\re{{\rm e}}        % number e
\def\rd{{\rm d}}        % differential

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\fC{{\mathfrak C}}
\def\fL{{\mathfrak L}}
\def\cL{{\mathcal L}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\bQ{{\mathbf Q}}
\def\bJ{{\mathbf J}}

\def\bV{{\mathbf V}}
\def\bA{{\bf A}}
\def\bB{{\bf B}}
\def\bC{{\bf C}}
\def\bc{{\mathbf c}}
\def\bb{{\mathbf b}}
\def\bN{{\mathbf N}}
\def\bR{{\mathbf R}}

\def\bGamma{{\mathbf \Gamma}}
\def\bPi{{\mathbf \Pi}}




\def\br{{\bf r}}
\def\x{\times}
\def\ox{\otimes}

\def\fA{{\mathfrak A}}
\def\fB{{\mathfrak B}}
\def\fF{{\mathfrak F}}
\def\fG{{\mathfrak G}}
\def\fS{{\mathfrak S}}

\def\mG{{\mathbb G}}
\def\cZ{{\mathcal Z}}
\def\fZ{{\mathfrak Z}}

\def\W{{\mathbb W}}
\def\mP{{\mathbb P}}

\def\om{{\ominus}}


\def\Var{{\bf Var}}

\def\bO{{\bf 0}}
\def\bone{{\bf 1}}
\def\bzero{{\bf 0}}

\def\sM{{\sf M}}




\def\bD{{\bf D}}

\def\bH{{\bf H}}
\def\bh{{\bf h}}

\def\bI{{\bf I}}
\def\bT{{\bf T}}
\def\bU{{\bf U}}
\def\bu{{\bf u}}
\def\bS{{\bf S}}
\def\cF{{\cal F}}
\def\cW{{\mathcal W}}
\def\fW{{\mathfrak W}}
\def\cX{{\mathcal X}}
\def\fX{{\mathfrak X}}

\def\mF{{\mathbb F}}
\def\cK{{\cal K}}


\def\cJ{{\cal J}}
\def\cD{{\mathcal D}}

\def\cM{{\cal M}}
\def\cV{{\cal V}}
\def\cC{{\cal C}}
\def\cR{{\cal R}}

\def\sE{{\sf E}}
\def\sP{{\sf P}}

\def\cG{{\mathcal G}}
\def\cI{{\cal I}}
\def\cP{{\cal P}}
\def\cQ{{\cal Q}}
\def\fY{{\frak Y}}
\def\cY{{\cal Y}}
\def\cA{{\cal A}}
\def\cB{{\cal B}}
\def\cE{{\mathcal E}}
\def\cov{{\bf cov}}
\def\var{{\bf var}}



\def\cN{{\cal N}}
\def\cS{{\mathcal S}}
\def\cT{{\cal T}}

\def\bG{{\mathbf G}}
\def\bL{{\mathbf L}}


\def\F{{\mathbb F}}
\def\mM{{\mathbb M}}
\def\mX{{\mathbb X}}
\def\mY{{\mathbb Y}}
\def\mV{{\mathbb V}}

\def\im{\mbox{${\rm im\,}$}}


\def\mS{{\mathbb S}}
\def\mT{{\mathbb T}}
\def\mZ{{\mathbb Z}}





\def\eps{\varepsilon}
\def\Ups{\Upsilon}
\def\phi{\varphi}


\def\diag{\mathop{\rm diag}}    % diagonal matrix
\def\blockdiag{\mathop{\rm blockdiag}}    % diagonal matrix









%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\input{epsf}


\begin{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\title{\Large\bf
Hardy-Schatten Norms of Systems, Output Energy
Cumulants and Linear Quadro-Quartic  Gaussian
    Control}
%==============================================================================
\author{
Igor G. Vladimirov,
\qquad Ian R. Petersen
\thanks{This work is supported by the Australian Research Council.
The authors are with the School of Engineering and Information
Technology, University of New South Wales at
    the  Australian  Defence Force Academy, Canberra ACT 2600,
    Australia. E-mail: 
        {\tt igor.g.vladimirov@gmail.com,
        i.r.petersen@gmail.com}.}
}
%==============================================================================
\onecolumn
\pagestyle{plain}

\maketitle
\thispagestyle{empty}
%\pagestyle{empty}

%==============================================================================
\begin{abstract}
This paper is concerned with linear stochastic control systems in state space. The integral of the squared norm of the system output over a bounded time interval is interpreted as energy. The cumulants of the output energy in the infinite-horizon limit are related to Schatten norms of the system in the Hardy space of transfer functions and the risk-sensitive performance index. We employ a novel performance criterion which seeks to minimize a combination of the average value and  the variance of the output energy of the system per unit time. The resulting  linear quadro-quartic Gaussian control problem involves the $\cH_2$ and $\cH_4$-norms of the closed-loop system. We obtain equations for the optimal controller and outline a homotopy method which reduces the solution of the problem to the numerical integration of  a differential equation initialized by the standard linear quadratic Gaussian  controller.
\end{abstract}
%==============================================================================



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\paragraph*{Keywords and phrases:}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%linear Gaussian control, output energy cumulants, Schatten norms,
%quadro-quartic performance index.
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\paragraph*{2010 MSC codes:}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%93C05,          % Linear systems
%%93C35,          % Multivariable systems
%93D25,          % Input-output approaches
%93E20,          % Optimal stochastic control
%%93B52,         % Feedback control
%60H10.          % Stochastic ordinary differential equations






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\tableofcontents
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section{Introduction}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%\begin{itemize}
%\item
%definition of $\|\cdot\|_4$
%
%\item
%connection with the variance of output energy
%
%\item
%connection with the risk-sensitive performance criterion
%
%\item
%higher order Hardy norms $\|\cdot\|_{2k}$ and their connection with
%cumulants of output energy
%
%
%\item
%invariance of $\|\cdot\|_{2k}$ under taking the dual
%
%
%\item
%state-space representations of the $\cH_4$-norm
%
%
%
%
%\item
%appendix: verification of the norm property (convexity)
%
%\item
%appendix: second moments of bilinear forms of Gaussian random
%vectors
%\end{itemize}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

This paper is concerned with linear multi-input multi-output control systems, governed in state space by Ito stochastic differential equations, driven by a standard Wiener process which is regarded  as a random disturbance. The integral of the squared Euclidean norm of the system output over a bounded time interval  is interpreted as \textit{energy}. In the disturbance attenuation paradigm, the output energy is to be minimized  in some sense.

Linear Quadratic Gaussian (LQG) control \cite{AM_1990}, for example, seeks to  minimize the expectation of the output energy which, in the infinite-horizon limit, reduces to the squared $\cH_2$-norm of the closed-loop system in an appropriate Hardy space of transfer functions. An alternative performance index is employed in the Risk-Sensitive  and Minimum Entropy control theories \cite{MG_1991}. They utilise the expected value of the exponential of the output energy multiplied by a scaling parameter to adjust the risk sensitivity.  Risk-sensitive control  extends the LQG approach and is robust with respect to  Kullback-Leibler relative entropy bounded uncertainties in the random noise \cite{DJP_2000}.

The risk-sensitive performance index can be represented as a series expansion with respect to the energy scaling parameter. The coefficients of this series are the rates of the asymptotically linear growth of the cumulants of the output energy in the infinite-horizon limit. The cumulant growth rates are directly related to higher-order Schatten norms \cite{Simon_2005} of the transfer function of the system in an appropriate Hardy space. This allows the risk-sensitive criterion to be  viewed  as a linear combination of powers of Hardy-Schatten norms of the system whose weights are governed by the risk-sensitivity parameter in a very specific way. The ``reverse engineering''  of the risk-sensitive index suggests a wide family of performance criteria in the form of linear combinations of powers of the Hardy-Schatten norms. This gives rise to a class of \textit{output energy cumulant} (OEC) control problems which extend the risk-sensitive paradigm. In fact, the LQG approach can be considered to explore this freedom to a certain degree by retaining the first term (the squared $\cH_2$-norm of the system)  of the risk-sensitive index expansion.

The present paper develops the OEC control idea,  outlined above, by employing a performance criterion which seeks to minimize a combination of the average value and  the variance of the output energy of the system per unit time. The resulting linear quadro-quartic Gaussian (LQQG) control problem utilizes a \textit{quadro-quartic functional} as a finer truncation of the risk-sensitive performance index which retains the $\cH_2$  and $\cH_4$-norms of the closed-loop system and the risk-sensitive parameter.

The $\cH_4$-norm, which involves the Schatten 4-norm of matrices \cite{HJ_2007} and is referred to as the \textit{quartic norm}, was introduced in \cite{VKS_1996}  as a subsidiary construct in the anisotropy-based robust control theory for discrete-time stochastic systems. In the present study, the quartic norm plays a central role and, in addition to providing the next  term in the risk-sensitive index expansion, quantifies (via the $\cH_4$ to $\cH_2$-norms ratio) the time scale beyond which the infinite-horizon LQG cost starts manifesting itself in sample paths of the output energy of the system.  %Minimizing the quartic norm at a given level of the $\cH_2$-norm of the system would therefore shorten the time scales needed for the infinite-horizon LQG cost to be a relevant measure of the disturbance attenuation capabilities of the system.


We consider the LQQG problem in the class of linear stabilizing controllers with the same state dimension as the underlying plant. This allows equations for an optimal controller to be obtained by using  Frechet derivatives of the quadro-quartic performance index of the closed-loop system with respect to the state-space realization matrices of the controller. The resulting set of equations depends on the risk sensitivity parameter and yields the standard LQG controller for a zero value of the parameter.
We outline a homotopy method which regards the parameter as a fictitious time variable and reduces the solution of the set of equations to a problem involving  the numerical integration of  an ordinary  differential equation (ODE) initialized by the standard LQG controller.


In addition to its possible extension to the
discrete-time case, the LQQG approach may also find application in the control of quantum
stochastic  systems as an alternative to the risk-sensitive control
paradigm.



%We consider an optimal control
%problem for a linear time invariant multi-input multi-output system
%with a transfer function $F(s)$, a standard Wiener process $W$ at
%the input, and a Gaussian random process $Z$ at the output.
%%We use a novel {\it quadro-quartic}
%%performance criterion which is obtained by truncating that from the
%%Risk-Sensitive Control.
%Following the disturbance attenuation paradigm,
%\begin{equation}
%\label{ET}
%    \cE_T
%    :=
%    \int_{0}^{T}
%    |z_t|^2
%    \rd t
%\end{equation}
%is interpreted as the output energy of the system over the time
%interval $[0,T]$, which is to be minimized in a stochastic sense.
%With $\theta$ a positive parameter, we employ a novel {\it
%quadro-quartic} performance criterion
%\begin{equation}
%\label{QQ}
%    {\bQ}_{\theta}(F)
%    :=
%    \|F\|_2^2 + \frac{\theta}{2} \|F\|_4^4
%\end{equation}
%which involves not only the average value of the output energy,
%quantified by the $\cH_2$-norm $\|F\|_2$ of the system as $
%    \lim_{T\to \infty}
%    (\bE \cE_T/T)
%    =
%    \|F\|_2^2
%    =
%    (2\pi)^{-1}
%    \int_{-\infty}^{\infty}
%    \Tr
%    S(\omega)
%    \rd \omega
%$ and used in the LQG control,  but also the variance of this energy
%per unit time:
%\begin{equation}
%\label{varET}
%    \lim_{T\to \infty}
%    \frac{\var(\cE_T)}{T}
%    =
%    2\|F\|_4^4
%    =
%    \frac{1}{\pi}
%    \int_{-\infty}^{\infty}
%    \Tr
%    (
%        S(\omega)^2
%    )
%    \rd \omega.
%\end{equation}
%Here, $
%    S(\omega)
%    :=
%    \wh{F}(\omega) \wh{F}(\omega)^*
%$ is the spectral density of $Z$, and $\wh{F}(\omega):= F(i\omega)$,
%where $F(s)$  is assumed to be analytic in the complex right
%half-plane and vanishing at infinity. The $\cH_4$-norm $\|F\|_4$ in
%(\ref{varET}) utilizes the Schatten 4-norm of matrices
%\cite{HJ_2007,Simon_2005} and was introduced in \cite{VKS_1996}  as
%a subsidiary construct in the anisotropy-based robust control theory
%for discrete-time systems. In the present paper, $\|F\|_4$ plays a
%central role and quantifies the rate of the mean square convergence
%$\cE_T/T \to \|F\|_2^2$, so that, in view of (\ref{varET}), the root
%mean square deviation $\sqrt{\var (\cE_T)}$ of the output energy
%becomes negligible compared to its mean $\bE \cE_T$ on time scales
%$T\gg
%    (
%        \|F\|_4/\|F\|_2
%    )^4
%$. Generalizing the $\cH_2$ and $\cH_4$-norms,
%\begin{equation}
%\label{H2k}
%    \|F\|_{2k}
%     :=
%    \sqrt[2k]{
%    \frac{1}{2\pi}
%    \int_{-\infty}^{\infty}
%    \Tr(S(\omega)^k)
%    \rd \omega}
%    =
%    \sqrt[2k]{
%    \frac{1}{(2k-2)!!}
%    \lim_{T\to \infty}
%    \frac{\bK_k(\cE_T)}{T}}
%\end{equation}
%defines a Hardy-Schatten $2k$-norm of the system $F$ for any
%positive integer $k$. The second equality in (\ref{H2k}) links this
%norm with the asymptotic behaviour of the $k$th cumulant
%$\bK_k(\cE_T)$ of the output energy (\ref{ET}). For small $\theta$,
%the {\it quartic} norm $\|F\|_4$ provides the next correction to the
%$\cH_2$-norm in the series expansion of the expected
%exponential-of-quadratic functional
%$$
%    \frac{2}{\theta}
%    \lim_{T\to \infty}
%        \frac{\ln \bE\re^{\theta\cE_T/2}}{T}
%     =
%    -\frac{1}{2\pi \theta}
%    \int_{-\infty}^{\infty}
%    \ln\det(I-\theta S(\omega))
%    \rd \omega
%    =
%    \frac{1}{\theta}
%    \sum_{k= 1}^{\infty}
%    \frac{\theta^k}{k}
%    \|F\|_{2k}^{2k}
%     =
%    \bQ_{\theta}(F)
%    +
%    O(\theta^2),
%$$
%which is used, with $\theta < \|F\|_{\infty}^{-2}$, in the
%risk-sensitive and minimum entropy control theories \cite{MG_1991}.
%We show that for a strictly proper system $F:=(A,B,C)$, with state
%$X$ governed by an Ito SDE as
%\begin{equation}
%\label{Fstate}
%    \rd x_t
%     =
%     Ax_t\rd t + B\rd w_t,
%    \quad
%    z_t
%     =
%    Cx_t,
%\end{equation}
%where $A$, $B$, $C$ are appropriately dimensioned real matrices and
%$A$ is Hurwitz, the quartic norm is computed as
%\begin{equation}
%\label{H4state}
%    \|F\|_4^4
%      =
%    2
%    \Tr
%    (
%        CPQPC^{\rT}
%    )
%    =
%    2
%    \Tr
%    (
%        B^{\rT}QPQB
%    )
%     =
%    -4
%    \Tr
%    (
%        A(PQ)^2
%    )
%\end{equation}
%via the controllability and observability Gramians $P$ and $Q$ of
%the system (\ref{Fstate}) satisfying the algebraic Lyapunov
%equations $AP + PA^{\rT} + BB^{\rT} = 0$ and $A^{\rT}Q + Q A  +
%C^{\rT}C = 0$, so that the quadro-quartic functional (\ref{QQ})
%takes the form
%\begin{equation}
%\label{QQstate}
%    \bQ_{\theta}(F)
%    =
%    -
%    2
%    \Tr
%    (
%        APQ
%        (I+\theta PQ)
%    ).
%\end{equation}
%We then consider a plant with the standard Wiener process $W$ as the
%input noise, control input $U$, controlled output $Z$, observation
%output $Y$ and $n$-dimensional state $X$, governed by
%\begin{equation}
%\label{plant}
%    \rd x_t
%     =
%    A x_t\rd t  + B_1 \rd w_t + B_2 u_t\rd t,
%\quad
%    z_t
%     =
%    C_1 x_t + D_{12} u_t,
%\quad
%    \rd y_t
%     =
%    C_2 x_t\rd t + D_{21} \rd w_t,
%\end{equation}
%where $A$, $B_k$, $C_j$, $D_{jk}$ are appropriately dimensioned real
%matrices. The control signal $U$ is produced by a controller
%$K:=(a,b,c)$ with $n$-dimensional state $\Xi$ as
%\begin{equation}
%\label{y_to_xi_yxi_to_u}
%    \rd \xi_t
%     =
%    a\xi_t\rd t + b\rd y_t,
%    \quad
%    u_t
%     =
%    c\xi_t.
%\end{equation}
%A linear quadro-quartic Gaussian (LQQG) control problem for the
%plant (\ref{plant}) is formulated as minimizing the quadro-quartic
%functional $\bQ_{\theta}(F)$ of the closed-loop system
%\begin{equation}
%\label{F}
%    F
%    :=
%    (\cA, \cB, \cC)
%    =
%    \left(
%        \left[
%            \begin{array}{cc}
%                      a &  bC_2\\
%                      B_2 c & A
%            \end{array}
%        \right],\,
%        \left[
%            \begin{array}{c}
%                      bD_{21}\\
%                      B_1
%            \end{array}
%        \right],\,
%        \left[
%            \begin{array}{cc}
%                      D_{12} c & C_1
%            \end{array}
%        \right]
%    \right)
%\end{equation}
%with the $2n$-dimensional combined state $(\Xi,X)$ over the
%controller matrices $a$, $b$, $c$  from (\ref{y_to_xi_yxi_to_u})
%which make  $\cA$ Hurwitz:
%\begin{equation}
%\label{LQQG}
%    \bQ_{\theta}(F)
%    =
%    -
%    2
%    \Tr
%    (
%        \cA\cP\cQ
%        (I+\theta \cP\cQ)
%    )
%    \longrightarrow
%    \min,
%    \quad
%    K := (a,b,c)\
%    {\rm stabilizes}\
%    F,
%\end{equation}
%where $\cP$ and $\cQ$ are the controllability and observability
%Gramians of the closed-loop system (\ref{F}).
%%The LQQG controller seeks a
%%compromise between minimizing the mean value and the variance of the
%%output energy of the system and reduces to the standard LQG
%%controller in the limiting case $\theta = 0$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Variance of output energy and quartic norm \label{sec:H4}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Suppose $W := (w_t)_{t\in \mR}$ is a $m$-dimensional standard Wiener
process (initialised in the infinitely distant past) at
the input of a linear time invariant (LTI) system $F$
with a square integrable $\mR^{p\x m}$-valued impulse response function $f:= (f_t)_{t\> 0}$; see Fig.~\ref{fig:zfw}.
%==============================================================================
\begin{figure}[htb]
\begin{center}
\unitlength=1.0mm
%\linethickness{0.6pt}
\begin{picture}(50.00,10.00)
    \put(20,0){\framebox(10,10)[cc]{$F$}}
    \put(40,5){\vector(-1,0){10}}
    \put(20,5){\vector(-1,0){10}}
    \put(42,5){\makebox(0,0)[lc]{$W$}}
    \put(8,5){\makebox(0,0)[rc]{$Z$}}
\end{picture}
\caption{
    An LTI system $F$ with input $W$ and
    output $Z$.}
\label{fig:zfw}
\end{center}%\vskip-5mm
\end{figure}
%==============================================================================
The output $Z := (z_t)_{t\in \mR}$ of the system is a $\mR^p$-valued Gaussian random process  defined by the Ito stochastic integral
$    z_t
    :=
    \int_{-\infty}^{t}
    f_{t-s} \rd w_s
$.
 The mean
value of $Z$ is zero and the  covariance function is
\begin{equation}
\label{cov}
    c_t
    :=
    \bE(z_t z_0^{\rT})
    =
    \int_{0}^{+\infty}
    f_{s+t}f_s^{\rT}
    \rd
    s
    =
    \frac{1}{2\pi}
    \int_{-\infty}^{+\infty}
    S(\omega)
    \re^{i\omega t}
    \rd \omega
    =
    c_{-t}^{\rT},
    \qquad
    t\> 0,
\end{equation}
where
\begin{equation}
\label{S}
    S(\omega)
    :=
    \wh{F}(\omega) \wh{F}(\omega)^*
    =
    \int_{-\infty}^{+\infty}
    c_t
    \re^{-i \omega t}\
    \rd t
\end{equation}
is the spectral density of $Z$. Here, $(\cdot)^* :=
(\overline{(\cdot)})^{\rT}$ denotes the complex conjugate transpose of a
matrix, and
$    \wh{F}(\omega)
    :=
    F(i\omega)
    =
    \int_{0}^{+\infty}
    f_t
    \re^{-i \omega t}
    \rd t
$
is the Fourier transform of the impulse response, that is, the
boundary value of the transfer function of the system
$    F(v)
    :=
    \int_{0}^{+\infty}
    f_t
    \re^{-vt}
    \rd t
$, with $
    \Re v >0
$.
With $f$ assumed to be square integrable, $F$ belongs to the Hardy space
$\cH_2^{p\x m}$ of $\mC^{p\x m}$-valued functions of a complex
variable, analytic in the right half-plane and endowed with the
$\cH_2$-norm
\begin{equation}
\label{H2}
    \|F\|_2
    :=
    \sqrt{
        \int_{0}^{+\infty}
        \|f_t\|^2
        \rd t
    }
%    =
%    \sqrt{
%        \Tr c_0
%    }
    =
    \sqrt{
        \frac{1}{2\pi}
        \int_{-\infty}^{+\infty}
        \|\wh{F}(\omega)\|^2
        \rd
    \omega}.
\end{equation}
Here, the Plancherel theorem is used, and  $\|M\|:= \sqrt{\bra M,M\ket}$ denotes the Frobenius norm of a matrix $M$ generated by the inner product $\bra M, N\ket:= \Tr (M^*N)$, so that
$\|\wh{F}(\omega)\|^2 = \Tr S(\omega)$ is the trace of the spectral
density from (\ref{S}). In view of (\ref{cov}), $\|F\|_2^2 = \Tr c_0 =
\bE(|z_t|^2)$ is the variance of the output signal for any $t$.
For a finite
time horizon $T>0$, the random variable
\begin{equation}
\label{ET}
    \cE_T
    :=
    \int_{0}^{T}
    |z_t|^2
    \rd t
\end{equation}
is interpreted as the {\it output energy} of the system $F$ over the
time interval $[0,T]$, and
\begin{equation}
\label{eT}
    \epsilon_T
    :=
    \cE_T/T
\end{equation}
is the corresponding output energy rate. The
mean value of $\epsilon_T$ coincides with the squared $\cH_2$-norm of the system (\ref{H2}):
$
    \bE \epsilon_T
    =
    \|F\|_2^2
$.
This ensemble average can manifest itself
in sample paths of $\epsilon_T$ only by virtue of the law
of large numbers, provided $T$ is large
enough. Under additional assumptions on the system
$F$, the rate of the mean square convergence
$    \mathop{\rm l.i.m.}_{T\to +\infty}
    \epsilon_T
    =
    \| F\|_2^2
$ is quantified by the asymptotic
behaviour of the variance of $\epsilon_T$. The convergence rate is described by the lemma below in terms of the quantity
\begin{equation}
\label{H4}
    \|F\|_4
     :=
    \sqrt[4]{
    \frac{1}{2\pi}
    \int_{-\infty}^{+\infty}
    \|S(\omega)\|^2
    \rd \omega}
%    =
%    \sqrt[4]{
%    \int_{-\infty}^{\infty}
%    \|c_t\|^2
%    \rd t
%    }
    =
    \sqrt[4]{
    2
    \int_{0}^{+\infty}
    \|c_t\|^2
    \rd t
    }.
\end{equation}
This is a
continuous-time counterpart of the $\cH_4$-norm introduced
 as a subsidiary construct in the anisotropy-based
robust control of discrete-time systems \cite{VKS_1996}. The second equality in (\ref{H4}) follows from the Plancherel theorem
applied to the spectral density (\ref{S}). The systems $F$ with
$\|F\|_4<+\infty$ form a normed space $\cH_4^{p\x m}$. The
integrand $\|S(\omega)\|^2 = \Tr ((\wh{F}(\omega) \wh{F}(\omega)^*)^2)$ in (\ref{H4}) is the
fourth power of the Schatten 4-norm \cite[p.~441]{HJ_2007} of the
matrix $\wh{F}(\omega)$; see also \cite{Simon_2005}. The
$\cH_4$-norm $\|F\|_4$ will be referred to as the \textit{quartic norm} of the system $F$.

%Note that (\ref{H4}) is slightly different from the straightforward
%extension of the conventional norm in the Hardy space $\cH_4$ of
%scalar-valued functions to the matrix case
%$$
%    \sqrt[4]{
%    \frac{1}{2\pi}
%    \int_{-\pi}^{\pi}
%    \|\wh{F}(\omega)\|^4
%    \rd \omega}.
%$$
%Although both definitions coincide in the case, where at least one
%of the dimensions $p$ or $m$ equals 1, and yield topologically
%equivalent norms in the general MIMO setting by the inequalities
%$$
%    \frac{(\Tr \Sigma)^2}{r}
%    \<
%    \|\Sigma\|^2
%    \<
%    (\Tr \Sigma)^2
%$$
%for any positive semi-definite Hermitian matrix $\Sigma$ of rank
%$r$, the definition (\ref{H4}) is preferable due to the
%probabilistic interpretation below.

%====================================================================
\begin{lemma}
\label{lem:H4} Let $F \in \cH_2^{p\x m}\bigcap \cH_4^{p\x m}$. Then the variance of the
output energy rate (\ref{eT}) of the system behaves asymptotically
as%\vskip-5mm
\begin{equation}
\label{vareTasy}
    \var(\epsilon_T)
    \sim
    2 \|F\|_4^4/T,
    \qquad
    T\to +\infty.
\end{equation}
\end{lemma}
%====================================================================

\begin{proof}
 By applying
Lemma~\ref{lem:bilinear} of Appendix~A to the Gaussian random
vectors $z_s$ and $z_t$ and using (\ref{cov}), it follows that
$
    \cov(|z_s|^2, |z_t|^2)
    =
    2
    \|c_{s-t}\|^2
$.
Hence, the variance of the output energy (\ref{ET})  can be computed as
\begin{equation}
\label{varET}
    \var(\cE_T)
     =
    \int_{[0,T]^2}
    \cov(|z_s|^2, |z_t|^2)
    \rd s\rd t
     =
    2
    \int_{[0,T]^2}
    \|c_{s-t}\|^2
    \rd s\rd t
    =
    4T
        \int_{0}^{T}
        (
            1-u/T
        )
        \|c_u\|^2
        \rd u,
%        \\
%
%    & \sim &
%    4T
%        \int_{0}^{+\infty}
%        \|c_u\|^2
%        \rd u,
%    \quad
%    T \to +\infty.
\end{equation}
where use is made of the property $c_t=c_{-t}^{\rT}$ and the
invariance of the Frobenius norm of a matrix under the transpose.
Since the assumption $F \in \cH_4^{p\x m}$ ensures the square
integrability of the covariance function (\ref{cov}), then
\begin{equation}
\label{varET1}
    \lim_{T\to +\infty}
        \int_{0}^{T}
        (
            1-u/T
        )
        \|c_u\|^2
        \rd u
        =
        \int_{0}^{+\infty}
        \|c_u\|^2
        \rd u
\end{equation}
holds by Lebesgue's dominated convergence
theorem. Since $
    \d_T
        \int_{0}^{T}
        (
            1-u/T
        )
        \|c_u\|^2
        \rd u
        =
        T^{-2}
        \int_{0}^{T}u\|c_u\|^2\rd u
        \> 0
$, the convergence is monotonic. Now,
(\ref{vareTasy}) is obtained by using (\ref{eT}) and combining
(\ref{varET}) and (\ref{varET1}) with (\ref{H4}):
$
    \var(\epsilon_T)
    =
    \var(\cE_T)/T^2
    \sim
    4
        \int_{0}^{+\infty}
        \|c_u\|^2
        \rd u/T
        =
        2\|F\|_4^4/T
$ as $T\to +\infty$.
\end{proof}
%==============================================================================


%By combining (\ref{outputenergyrate}) and (\ref{limeT}) with
%(\ref{vareTasy}), it follows that
%\begin{equation}
%\label{varbarET}
%    \bE((\epsilon_T - \|F\|_2^2)^2)
%    =
%    \var(\epsilon_T)
%    \sim
%    \frac{2\|F\|_4^4}{T},
%    \quad
%    T \to +\infty.
%\end{equation}

In view of a central limit theorem for quadratic functionals of Gaussian processes
\cite[Theorem~2]{Ginovian_1994}, the relation (\ref{vareTasy}) provides the scaling factor for the asymptotic
standard normality of the random variable $\sqrt{T/2} (\epsilon_T - \|F\|_2^2)/\|F\|_4^2$ as $T\to +\infty$. Heuristically,  the root mean square deviation of $\epsilon_T$ from its mean
value $\|F\|_2^2$ is relatively small if
\begin{equation}
\label{T*}
    T
    \gg
    T_*
    :=
    2
    (
        \|F\|_4/\|F\|_2
    )^4.
\end{equation}
The right-hand side of (\ref{T*}) quantifies the time horizon beyond
which the $\cH_2$-norm $\|F\|_2$ manifests itself in the sample
paths of the output energy of the system. On the other hand, for $T
\ll T_*$, the ergodic properties of the system output $Z$ do not
expose themselves since the expected value $\bE \epsilon_T =
\|F\|_2^2$ of the output energy rate is ``indistinguishable'' in the
background of random fluctuations whose standard deviation can be estimated by using
(\ref{vareTasy}) as
$
    \sqrt{\var(\epsilon_T)}
    \sim
    \|F\|_4^2\,
    \sqrt{2/T}
    \gg
    \bE \epsilon_T
$. Thus, the squared $\cH_2$-norm as the
average output energy loses its significance for quantifying the disturbance attenuation capabilities of the system on short time scales $T\ll T_*$. The critical time horizon $T_*$ defined by (\ref{T*})
is similar to the {\it integral time scale} of measurements  in
turbulent flows \cite[pp.~50--51]{Frisch_1995}. As an example, let $Z$ be an Ornstein-Uhlenbeck
process generated from a standard Wiener process $W$ by a single-input single-output system $F$ according to the SDE
\begin{equation}
\label{OU}
    \rd z_t = az_t\rd t + \sqrt{2|a|}\,\rd w_t,
\end{equation}
parameterized by $a <0$.
The covariance function (\ref{cov}) of $Z$ is
$
    c_t
    =
    \re^{a|t|}
$, and
the $\cH_2$
and $\cH_4$-norms of the system $F$, defined by (\ref{H2}) and (\ref{H4}),  are
$
    \|F\|_2 = 1
$
and
$
    \|F\|_4
    =
    |a|^{-1/4}
$.
Therefore, the critical time
horizon (\ref{T*}) takes the form  $    T_*
    =
    2/|a|
$  and coincides with the typical transient time of the process; see Fig.~\ref{fig:energy}.
%==============================================================================
\begin{figure}[htb]
\begin{center}
\includegraphics[width=80mm]{energy.ps}
\caption{
    100 sample paths of $\epsilon_T$ versus $T\<10$ for the Ornstein-Uhlenbeck process
    generated by  (\ref{OU}) with $a=-1$, so that the critical time horizon
     beyond which $\epsilon_T$ exposes relative proximity
    to the limit value $\|F\|_2^2 = 1$ (horizontal bold line) is $T_*= 2$.
    The dashed bold lines localize the typical values of $\epsilon_T$ which form a ``tube'' of half-width
    $\sqrt{T_*/T}$ about the limit.
}
    \label{fig:energy}
\end{center}
\end{figure}
%==============================================================================


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Cumulants of output energy and Hardy-Schatten norms\label{sec:cumulants}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

For a finite time horizon $T>0$, let $C_T$ denote a Toeplitz  integral
operator whose kernel is specified by the covariance function
(\ref{cov}). An $\mR^p$-valued integrable function $\psi:=
(\psi_t)_{0\<t \<T}$ is mapped by $C_T$ to $\phi:= (\phi_s)_{0\< s\< T}$ as
$    \phi_s
    :=
    \int_{0}^{T}
    c_{s-t}\psi_t
    \rd
    t
$.
Suppose $\theta$ is a real parameter
satisfying $0<\theta < 1/\rho(C_T)$, where $\rho(\cdot)$ is the spectral radius. In view of the Fredholm
formula \cite[Theorem~3.10 on p.~36]{Simon_2005} (see also \cite{Ginovian_1994} and references therein),
\begin{equation}
\label{expquad}
    \ln \bE
    \re^{\theta\cE_T/2}
     =
    -
    \frac{1}{2}
    \Tr
    \ln(
        I-\theta C_T
    )
     =
        \frac{1}{2}
        \sum_{k\> 1}
        \theta^k
        \Tr(C_T^k)/k,
\end{equation}
where $I$ is the identity operator. The trace of the $k$-fold iterate
of $C_T$ is computed as
\begin{equation}
\label{TrCT}
    \Tr(C_T^k)
    =
    \int_{[0,T]^k}
    \Tr
    (
        c_{t_0-t_1}
        c_{t_1-t_2}
        \x
        \ldots
        \x
        c_{t_{k-2}-t_{k-1}}
        c_{t_{k-1}-t_0}
    )
    \rd t_0
    \x \ldots \x
    \rd t_{k-1}.
\end{equation}
The expectation in (\ref{expquad}) is the
moment-generating function of $\cE_T$, and
hence,
\begin{equation}
\label{expquad1}
    \ln \bE\re^{\theta\cE_T/2}
    =
    \sum_{k\> 1}
    (\theta/2)^k
    \bK_k(\cE_T)/k!
    =
    \theta
    \left(
        \bE\cE_T
        +
        \theta
        \var(\cE_T)/4
    \right)/2
    +
    O(\theta^3),
    \qquad
    \theta \to 0.
\end{equation}
Here,
$
    \bK_k(\xi)
    :=
    \left.
        \d_v^k
        \ln \bE \re^{v\xi}
    \right|_{v=0}
    =
    P_k(\bE \xi, \ldots, \bE (\xi^k))
$
denotes the $k$th cumulant of a random variable $\xi$, which is
related with the first $k$ moments of $\xi$ via a universal
polynomial $P_k$. The first three of these
polynomials are
$
P_1(\mu_1)
=
    \mu_1
$,
$
P_2(\mu_1, \mu_2)
 =
    \mu_2 - \mu_1^2
$ and
$
P_3(\mu_1, \mu_2, \mu_3)
  =
    \mu_3 - 3\mu_1\mu_2 + 2\mu_1^3$.
By comparing the power series in (\ref{expquad}) and
(\ref{expquad1}) and using the identity $(2r)!! = r!2^r$, it follows that
the $k$th cumulant of the output energy (\ref{ET})
of the system is related to the trace (\ref{TrCT}) as
\begin{equation}
\label{Kzeta}
    \bK_k(\cE_T)
    =
    (2k-2)!!
    \Tr(C_T^k).
\end{equation}
Using (\ref{S}) and extending (\ref{H2}) and (\ref{H4}),
we define, for a positive  integer $k$, a
higher order Hardy norm of the system $F$ by
\begin{equation}
\label{H2k}
    \|F\|_{2k}
     :=
    \sqrt[2k]{
    \frac{1}{2\pi}
    \int_{-\infty}^{+\infty}
    \Tr(S(\omega)^k)
    \rd \omega},
\end{equation}
which reproduces the $\cH_2$
and $\cH_4$-norms for $k=1,2$. Here, $\sqrt[2k]{\Tr(S(\omega)^k)}$
is the Schatten $2k$-norm \cite[p.~441]{HJ_2007} of the matrix
$\wh{F}(\omega)$. The resulting  Hardy-Schatten
space $\cH_{2k}^{p\x m}$ is equipped with the
norm $\|\cdot\|_{2k}$.
%==============================================================================
Similarly to the $\cH_2$-norm, the $\cH_{2k}$-norms (\ref{H2k}) are all
invariant under replacing the system $F$ with its dual $F^{\dagger}$,
\begin{equation}
\label{dual_invariance}
    \|F^{\dagger}\|_{2k}
    =
    \|F\|_{2k},
    \qquad
    k \> 1,
\end{equation}
where $F^{\dagger}$ has the transposed impulse response
$(f_t^{\rT})_{t \> 0}$. Indeed, the transpose of a square matrix
does not modify its spectrum, and for conformable complex
matrices $X$ and $Y$, the matrices $XY$ and
$YX$ share nonzero eigenvalues. Therefore, with the
dependence  on the frequency $\omega$ omitted for brevity,
$
    \Tr
    (
        (
            \wh{F}^{\rT}
            (
                \wh{F}^{\rT}
            )^*
        )^k
    )
    =
    \Tr
    (
        (
            (
                \wh{F}^*
                \wh{F}
            )^k
        )^{\rT}
    )
    =
    \Tr
    (
        (
            \wh{F}
            \wh{F}^*
        )^k
    )
$,
and hence (\ref{dual_invariance}) follows.
%, so that they form a decreasing sequence
%$$
%    \cH_{\infty}^{p\x m}
%    \subset
%    \cH_{2k+2}^{p\x m}
%    \subset
%    \cH_{2k}^{p\x m} \subset \cH_2^{p\x m},
%    \quad
%    k> 1.
%$$
%Moreover, the norms (\ref{H2k}) are convergent  to the
%$\cH_{\infty}$-norm in the sense that for any $F \in
%\cH_{\infty}^{p\x m}$,
%$$
%    \lim_{k \to +\infty}
%    \|F\|_{2k}
%    =
%    \|F\|_{\infty}
%    :=
%    \mathop{\rm ess\, sup}_{-\pi \< \omega \< \pi}
%    \sigma_{\rm max}(\wh{F}(\omega)),
%$$
%where $\sigma_{\rm max}$ is the largest singular value of a matrix.
By the Szeg\H{o} limit theorem for Toeplitz operators
\cite{GS_1958}, under additional integrability conditions,
\begin{equation}
\label{Szego}
    \lim_{T \to +\infty}
    \frac{\Tr \chi(C_T)}{T}
    =
    \frac{1}{2\pi}
    \int_{-\infty}^{+\infty}
    \Tr \chi(S(\omega)) \rd \omega.
\end{equation}
Here, $\chi$ is a function of
a complex variable, satisfying $\chi(0) = 0$ and analytic in a neighbourhood of the
interval $[0,\|F\|_{\infty}^2]$, with $\|F\|_{\infty}$ the $\cH_{\infty}$-norm of $F$.
In view of (\ref{Kzeta}), the application of (\ref{Szego}) to elementary polynomials $\chi(v):=v^k$ yields
the asymptotically linear growth of the output energy cumulants with respect to time:
$    \lim_{T\to +\infty}
    (\bK_k(\cE_T)/T)
    =
    (2k-2)!!
    \|F\|_{2k}^{2k}
$,
provided $F\in \bigcap_{j=1}^{k} \cH_{2j}^{p\x m}$, with
Lemma~\ref{lem:H4} being a particular case for $k=2$.  The application of  (\ref{Szego})
to $\chi(v):= (2/\theta)\ln(1-\theta v)$, with
$0< \theta< \|F\|_{\infty}^{-2}$, gives
\begin{align}
\nonumber
    \frac{2}{\theta}\lim_{T\to +\infty}
        \frac{\ln \bE\re^{\theta\cE_T/2}}{T}
     & =
    -\frac{1}{2\pi\theta}
    \int_{-\infty}^{+\infty}
    \ln\det(I_p-\theta S(\omega))
    \rd \omega\\
\label{quartic_correction}
    & =
    \sum_{k\> 1}
    \theta^{k-1}
    \|F\|_{2k}^{2k}/k
     =
    \bQ_{\theta}(F)
    +
    O(\theta^2),
    \qquad
    \theta \to 0+,
\end{align}
where $I_p$ denotes  the identity matrix of order $p$, and
\begin{equation}
\label{QQ}
    \bQ_{\theta}(F)
    :=
    \|F\|_2^2
    +
    \theta\|F\|_4^4\big/2.
\end{equation}
The expected exponential-of-quadratic functional $\bE\re^{\theta\cE_T/2}$ in (\ref{quartic_correction}) is used as a
performance criterion in the risk-sensitive  and minimum entropy
control theories \cite{MG_1991}. The \textit{quartic norm} $\|F\|_4$
provides the next correction to the squared $\cH_2$-norm $\|F\|_2^2$ in the
series expansion (\ref{quartic_correction})  for small $\theta$.
Therefore, the {\it quadro-quartic} functional $\bQ_{\theta}$, defined by
(\ref{QQ}), can be regarded as a finer truncation of the risk-sensitive performance
index.






%==============================================================================
%We have also used the
%property that the Frobenius norm of a Hermitian matrix $\Sigma$
%satisfies $\|\Sigma\|^2 = \Tr(\Sigma^2)$, where the latter quantity
%is the sum of squared eigenvalues of $\Sigma$.


%$$
%    E_k = \frac{\rd^k}{\rd \theta^k} \bE\re^{\theta|\xi|^2}
%$$


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Quadro-quartic functional in state space}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Let $F$ be a strictly proper LTI system with an $m$-dimensional standard Wiener process
$W$ at the input, $p$-dimensional
output $Z$ and $n$-dimensional state $X$ governed by an Ito SDE:
\begin{equation}
\label{Fstate}
    \rd x_t
     =
     Ax_t\rd t + B\rd w_t,
    \qquad
    z_t
     =
    Cx_t,
\end{equation}
where $A\in \mR^{n\x  n}$,
$B\in \mR^{n\x m}$, $C\in \mR^{p\x n}$ are constant matrices. The state-space
representation  will be written as
%\vskip-5mm
\begin{equation}
\label{FABC}
    F
    =
    (A,B,C)
    =
        \begin{array}{rl}
         &   \!\!\!\!\!\!{}_{\leftarrow n \rightarrow}{}_{\leftarrow m \rightarrow} \\
            \begin{array}{cc}
                {}^n & \hskip-2mm\updownarrow \\
                {}_p & \hskip-2mm\updownarrow
            \end{array} &\hskip-4mm
            \left[
                \begin{array}{c|c}
                    A & B \\
                    \hline
                    C & 0 \\
                \end{array}
            \right]\\
            {}
        \end{array},
\end{equation}
where we have also shown the dimensions, and the horizontal and vertical separators serve to avoid
confusion with an ordinary block matrix.
%; see Fig.~\ref{fig:abcd}.
%%==============================================================================
%\begin{figure}[htb]
%\begin{center}
%\includegraphics[width=6cm]{Images/abcd.ps}
%\end{center}
%\caption{
%    \label{fig:abcd}
%    The system (\ref{Fstate}) with input $W$,
%    output $Z$ and state $X$.
%}
%\end{figure}
%%==============================================================================
The dual system is $
    F^{\dagger}
    =
    (A^{\rT}, C^{\rT},B^{\rT})
$. If the matrix $A$ is Hurwitz, then the mutually dual
controllability and observability Gramians $P$ and $Q$  of (\ref{FABC}) are
unique solutions of the algebraic Lyapunov
 equations
\begin{equation}
\label{PQ}
    AP + PA^{\rT} + BB^{\rT} = 0,
    \qquad
    A^{\rT}Q + Q A  + C^{\rT}C = 0.
\end{equation}
In what follows, an important role is played by the matrix
\begin{equation}
\label{H}
    H:= QP,
\end{equation}
whose spectrum is formed by the squared Hankel singular values of the system (\ref{FABC}).
We will write $\|X\|_M:= \sqrt{\Tr(X^{\rT} M X)}$ for the weighted Frobenius (semi-) norm of a
real matrix $X$ generated by a positive (semi-) definite matrix $M$.






%==============================================================================
\begin{lemma}
\label{lem:H4state} Let $F$ be an asymptotically stable system with
the  state-space realization (\ref{FABC}).  Then the quartic norm (\ref{H4})
 is expressed in terms of the Gramians $P$, $Q$ from (\ref{PQ}) and the matrix $H$ from (\ref{H}) as
\begin{equation}
\label{H4state}
    \|F\|_4^4
     =
    2
    \|(A,PC^{\rT},C)\|_2^2
%    =
%    \left\|
%        \left[
%            \begin{array}{c|c}
%            A & PC^{\rT}\\
%            \hline
%            C & 0
%            \end{array}
%        \right]
%    \right\|_2^2
    =
    2\|PC^{\rT}\|_Q^2
    =
    2
    \|(A,B,B^{\rT}Q)\|_2^2
%    \left\|
%        \left[
%            \begin{array}{c|c}
%            A & B\\
%            \hline
%            B^{\rT}Q& 0
%            \end{array}
%        \right]
%    \right\|_2^2
    =
    2\|QB\|_P^2
      =
    -4\Tr(A^{\rT}H^2).
\end{equation}
\end{lemma}
%==============================================================================

\begin{proof}
Let $Z$ be
a stationary Gaussian random process generated by
(\ref{Fstate}), with $W$ a standard Wiener process.  Then the
steady-state covariance function (\ref{cov}) is
\begin{equation}
\label{ct}
    c_t
    =
    C\re^{At} P C^{\rT},
    \qquad
    t\> 0.
\end{equation}
Here, we use the fact that the controllability
Gramian is the steady-state covariance matrix of the state of the system:
$P=\cov(x_t)$. Since the function
$c_t$ in (\ref{ct}) coincides with the impulse response
of the system $(A,PC^{\rT}, C)$, then (\ref{H4}) yields
$    \|F\|_4^4
     =
    2
    \|(A,PC^{\rT},C)\|_2^2
%    \left\|
%        \left[
%            \begin{array}{c|c}
%            A & PC^{\rT}\\
%            \hline
%            C & 0
%            \end{array}
%        \right]
%    \right\|_2^2
    =
    2
    \Tr
    (
        CPQPC^{\rT}
    )
    =
    2\|PC^{\rT}\|_Q^2
$,
which  proves the first two equalities in (\ref{H4state}). Here,
we have also used the property that the system $(A,PC^{\rT},C)$  shares the matrices $A$, $C$ with the underlying
system (\ref{FABC}) and hence, inherits from $F$ the observability
Gramian $Q$. The remaining three equalities in
(\ref{H4state}) follow from the first two by the invariance of the
$\cH_2$ and $\cH_4$-norms under taking the dual of a system, and by  the duality
of the controllability and observability Gramians.
\end{proof}

%%==============================================================================
%\textit{Remark~2:}
%It may seem counterintuitive that the linear equations (\ref{PQ})
%imply a substantially nonlinear identity
%\begin{equation}
%\label{comm}
%    \Tr(CPQPC^{\rT})
%    =
%    \Tr(B^{\rT} QPQB)
%\end{equation}
%which was established in the proof of Lemma~\ref{lem:H4state} using
%an indirect duality argument. Nevertheless, (\ref{comm}) is amenable
%to a direct verification. Right multiplication of both sides of the
%first equality from (\ref{PQ}) by $Q$ yields
%\begin{equation}
%\label{comm1}
%    APQ +
%    P\underbrace{A^{\rT}Q} +
%    BB^{\rT}Q
%    =
%    0.
%\end{equation}
%The Lyapunov equations (\ref{PQ}) can be regarded as relations which
%show how the Gramians commute (or, rather, anticommute) with the
%matrices $A$ and $A^{\rT}$. Indeed, the second equation in
%(\ref{PQ}) implies
%\begin{equation}
%\label{comm2}
%    A^{\rT}Q
%    =
%    -QA - C^{\rT}C
%\end{equation}
%which (up to the transpose of $A$) describes an anticommutation
%relation for the observability Gramian $Q$ with respect to $A$. This
%``anticommutation'' form of the Lyapunov equation allows the
%matrices $A$ or $A^{\rT}$ to be ``pushed'' through the product of
%alternating $P$ and $Q$. Substituting (\ref{comm2}) into
%(\ref{comm1}) yields
%\begin{equation}
%\label{comm3}
%    APQ - PQA - PC^{\rT}C + BB^{\rT}Q = 0.
%\end{equation}
%Therefore,  by left multiplication of both sides of (\ref{comm3}) by
%$PQ$ and taking the trace, it follows that
%\begin{multline}
%\label{comm3}
%    \Tr
%    (
%        PQAPQ - PQPQA - PQPC^{\rT}C + PQBB^{\rT}Q
%    )\\
%     =
%    -
%    \Tr(PQPC^{\rT}C)
%    +
%    \Tr(PQBB^{\rT}Q)
%    = 0.
%\end{multline}
%Here, the first equality follows from the cyclic property of the
%trace that $\Tr(MN) = \Tr(NM)$ whereby  $\Tr(PQAPQ) = \Tr(PQPQA)$.
%Application of the cyclic property to the last two traces in
%(\ref{comm3}) gives (\ref{comm}). \endremark
%
%
%
%%==============================================================================

The controllability and observability Gramians $\Phi$, $\Psi$ of a subsidiary system $(A,PC^{\rT}, B^{\rT}Q)$, which satisfy the algebraic Lyapunov equations
\begin{equation}
\label{Schattenian_Phi_Psi}
    A\Phi + \Phi A^{\rT} +  P C^{\rT}C  P
      =
    0,\qquad
    A^{\rT} \Psi
    +
    \Psi A +  Q BB^{\rT}  Q
      =  0,
\end{equation}
 will be referred to as the controllability and observability \textit{Schattenians} of the system (\ref{FABC}). The representations (\ref{H4state}) imply that
% \begin{equation}
% \label{Schattenian_quartic}
$$    \|F\|_4^4
    =
    2\Tr(C \Phi C^{\rT})
    =
    2\Tr(B^{\rT} \Psi B),
$$
% \end{equation}
and hence, the significance of the Schattenians $\Phi$, $\Psi$ for the quartic norm is analogous to the role which the Gramians $P$, $Q$ play for the $\cH_2$-norm.
%The following theorem is obtained by combining Lemma~\ref{lem:H4state} with the state-space representation of the $\cH_2$-norm.

%==============================================================================
\begin{theorem}
\label{th:QQstate} Let $F$ be an asymptotically stable system with
the  state-space realization (\ref{FABC}).  Then the quadro-quartic
functional (\ref{QQ}) is expressed in terms of the Gramians  $P$, $Q$  from (\ref{PQ}) and the matrix $H$ from (\ref{H}) as
\begin{align}
\nonumber
    \bQ_{\theta}(F)
    & =
    \left\|
        (
            A,
            \begin{bmatrix}
            B & \sqrt{\theta} PC^{\rT}
            \end{bmatrix},
            C
        )
    \right\|_2^2\\
\nonumber
%    \left\|
%        \left[
%            \begin{array}{c|cc}
%            A & B & \sqrt{\theta} PC^{\rT}\\
%            \hline
%            C & 0 & 0
%            \end{array}
%        \right]
%    \right\|_2^2\\
    & =
    \Tr
    (
        (
            BB^{\rT} + \theta PC^{\rT}CP
        )
        Q
    )\\
\nonumber
    & =
    \left\|
        \left(
            A,
            B,
            \begin{bmatrix}
                C \\
\nonumber
                 \sqrt{\theta} B^{\rT}Q
            \end{bmatrix}
        \right)
    \right\|_2^2\\
%    \left\|
%        \left[
%            \begin{array}{c|c}
%            A & B\\
%            \hline
%            C & 0\\
%            \sqrt{\theta} B^{\rT}Q & 0
%            \end{array}
%        \right]
%    \right\|_2^2\\
    \label{QQstate}
    & =
    \Tr
    (
        (
            C^{\rT}C + \theta QBB^{\rT}Q
        )
        P
    )
    =
    -
    2\Tr(A^{\rT}H(I_n+\theta H)).
\end{align}
\end{theorem}
%==============================================================================
\begin{proof}
Substitution of
$\|F\|_2 = \sqrt{\Tr (B^{\rT}QB)}$ and the first two equalities from
(\ref{H4state}) into (\ref{QQ}) yields
\begin{align*}
    \bQ_{\theta}(F)
    & =
    \|(A,B,C)\|_2^2
    +
    \theta
    \|(A,PC^{\rT},C)\|_2^2\\
    & =
    \|
        (
            A,
                \begin{bmatrix}
                    B & \sqrt{\theta}PC^{\rT}
                \end{bmatrix},
            C
        )
    \|_2^2
%    =
%    \Tr
%    (
%        BB^{\rT}Q
%    )
%    +
%    \theta
%    \Tr
%    (
%        CPQPC^{\rT}
%    )\\
    =
    \Tr
    (
        (BB^{\rT} + \theta PC^{\rT}CP)
        Q
    ),
\end{align*}
which establishes the first two equalities in (\ref{QQstate}).
 The third and fourth representations of the quadro-quartic
functional are obtained from the first two by the duality argument
or directly from the third and fourth equalities in (\ref{H4state}).
The last representation of $\bQ_{\theta}(F)$ in
(\ref{QQstate}) follows from the previous ones by using the Lyapunov
equations (\ref{PQ}):
\begin{align*}
    \bQ_{\theta}(F)
    & =
    \Tr
    (
        (BB^{\rT} + \theta PC^{\rT}CP)
        Q
    )\\
    & =
    -
    \Tr
    (
        (AP+PA^{\rT} + \theta P(A^{\rT}Q + QA)P)
        Q
    )\\
    & =
    -2
    \Tr
    (
        A^{\rT}QP
        +
        \theta
        A^{\rT}(QP)^2
    )
    =
    -2
    \Tr
    (
        A^{\rT}H
        (
            I_n
            +
            \theta
            H
        )
    ).
\end{align*}
\end{proof}
%==============================================================================

%\textit{Remark~3:} The arithmetic mean of the traces from (\ref{QQstate})
%provides a ``symmetrized'' representation  of the  quadro-quartic
%functional:
%\begin{equation}
%\label{symQQstate}
%    \bQ_{\theta}(F)
%    =
%    \frac{1}{2}
%    \Tr
%    (
%        C^{\rT}CP
%        +
%        BB^{\rT}Q
%        +
%        \theta
%        (
%            C^{\rT}CPQP
%            +
%            BB^{\rT}QPQ
%        )
%    ).
%\end{equation}
%





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Linear quadro-quartic Gaussian control problem}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Consider a plant with an  $m_1$-dimensional standard Wiener
process $W$ as the input disturbance and an $m_2$-dimensional input
control signal $U$. The outputs of the system are a
$p_1$-dimensional to-be-controlled signal $Z$ and a
$p_2$-dimensional observation signal $Y$.
Also, the system has an $n$-dimensional state $X$. These processes are governed by
\begin{eqnarray}
\label{nextx}
    \rd x_t
    & = &
    A x_t\rd t \ + \ B_1 \rd w_t + B_2 u_t\rd t,\\
\label{z}
    z_t
    & = &
    C_1 x_t \ \ \, + \qquad \quad\quad \ D_{12} u_t,\\
\label{y}
    \rd y_t
    & = &
    C_2 x_t\rd t + D_{21} \rd w_t.
\end{eqnarray}
Here,
$
    A\in \mR^{n\x n}$,
$
    B_k\in \mR^{n\x m_k}$,
$
    C_j\in \mR^{p_j\x n}
$,
$
    D_{jk}\in \mR^{p_j\x m_k}
$, with $D_{11}=0$ and
$D_{22}=0$. The control signal $U$ is generated at the output of a
controller $K$ with input $Y$. We consider a strictly proper LTI
controller
\begin{equation}
\label{K}
    K
    =
        \begin{array}{rl}
         &   \!\!\!\!\!{}_{\leftarrow n \rightarrow}{}_{\leftarrow p_2 \rightarrow} \\
            \begin{array}{cc}
                {}^n & \hskip-3mm\updownarrow \\
                {}_{m_2} & \hskip-3mm\updownarrow
            \end{array} &\hskip-4mm
            \left[
                \begin{array}{c|c}
                    a \ & \ b \, \\
                    \hline
                    c \ & \ 0 \,\\
                \end{array}
            \right]\\
            {}
        \end{array},
\end{equation}
with an $n$-dimensional state $\Xi$. It is driven by the observation
$Y$ and produces the output $U$ as
\begin{equation}
\label{y_to_xi_yxi_to_u}
    \rd \xi_t
     =
    a\xi_t\rd t + b\rd y_t,
    \qquad
    u_t
     =
    c\xi_t,
\end{equation}
where $a\in \mR^{n\x n}$, $b\in \mR^{n\x p_2}$, $c\in \mR^{m_2\x
n}$. The
closed-loop system
%\vskip-7mm
\begin{equation}
\label{F}
    F
    :=
        \begin{array}{rl}
         &   \!\!\!\!\!{}_{\leftarrow 2n \rightarrow}{}_{\,\leftarrow m_1 \rightarrow} \\
            \begin{array}{cc}
                {}^{2n} & \hskip-2mm\updownarrow \\
                {}_{p_1} & \hskip-2mm\updownarrow
            \end{array} &\hskip-4mm
            \left[
                \begin{array}{c|c}
                    \ \cA \ & \ \cB \, \\
                    \hline
                    \ \cC\ & \ 0\, \\
                \end{array}
            \right]\\
            {}
        \end{array}    
%    \left[
%        \begin{array}{c|c}
%            \cA & \cB\\
%            \hline
%            \cC & 0
%        \end{array}
%    \right]
    =
    \left[
        \begin{array}{cc|c}
                  a &  bC_2 & bD_{21}\\
                  B_2 c & A & B_1\\
            \hline
                  D_{12} c & C_1 & 0
        \end{array}
    \right],
\end{equation}
governed
by (\ref{nextx})--(\ref{y_to_xi_yxi_to_u}) and depicted in Fig.~\ref{fig:plant},
%==============================================================================
\begin{figure}[htb]
\begin{center}
\unitlength=1.0mm
\begin{picture}(50.00,25.00)
    \put(20,15){\framebox(10,10)[cc]{{\small plant}}}
    \put(40,23){\vector(-1,0){10}}
    \put(20,23){\vector(-1,0){10}}
    \put(42,23){\makebox(0,0)[lc]{$W$}}
    \put(8,23){\makebox(0,0)[rc]{$Z$}}
    \put(40,17){\vector(-1,0){10}}
    \put(20,17){\line(-1,0){10}}
    \put(10,17){\line(0,-1){12}}
    \put(10,5){\vector(1,0){10}}
    \put(30,5){\line(1,0){10}}
    \put(40,5){\line(0,1){12}}
    \put(42,11){\makebox(0,0)[lc]{$U$}}
    \put(8,11){\makebox(0,0)[rc]{$Y$}}
    \put(20,0){\framebox(10,10)[cc]{$K$}}
\end{picture}
\caption{
    The closed-loop system $F$ with input $W$ and
    output $Z$.}
\label{fig:plant}
\end{center}%\vskip-5mm
\end{figure}
%==============================================================================
has the $2n$-dimensional
combined state $(\Xi,X)$. We formulate a {\it linear
quadro-quartic Gaussian} (LQQG) control problem as the minimization of the
functional (\ref{QQ}) over $n$-dimensional controllers (\ref{K})
such that the matrix $\cA$ of the closed-loop system in
(\ref{F}) is Hurwitz:
\begin{equation}
\label{LQQG}
    \bQ:=\bQ_{\theta}(F)
    =
    -2\Tr(\cA^{\rT}H(I_{2n}+\theta H))
    \longrightarrow
    \min,
    \qquad
    K\
    {\rm stabilizes}\
    F.
\end{equation}
Here, $\theta\>0$ is a given parameter as before, and use is made of Theorem~\ref{th:QQstate}, so that the matrix $H$ is associated by (\ref{H}) with the Gramians $P$, $Q$ of the closed-loop system satisfying  the algebraic Lyapunov equations
\begin{equation}
\label{cPQ}
    \cA P
    +
     P \cA^{\rT}
    +
    \cB\cB^{\rT}
    = 0,
\qquad
    \cA^{\rT} Q
    +
     Q  \cA
    +
    \cC^{\rT}\cC
    = 0.
\end{equation}
In the case $\theta = 0$, the LQQG problem (\ref{LQQG}) reduces to the standard linear quadratic Gaussian (LQG) control problem. For $\theta > 0$, the LQQG problem is a compromise between minimizing the mean value and the variance of the output energy per unit time, with $\theta$ becoming  the relative weight of the quartic norm.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Matrices with $\Gamma$-shaped sparsity\label{sec:bGamma}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Since it is convenient to assemble the state-space realization matrices into a matrix with
``$\Gamma$-shaped'' sparsity,  we denote the set of real $(r+p)\x(r+m)$-matrices
with zero bottom-right block of size $(p\x m)$ by
%\vskip-3mm
\begin{equation}
\label{bGamma}
    \bGamma_{r,m,p}
    :=
    \left\{
                 \begin{bmatrix}
                     \rho  &  \sigma \\
                     \tau & 0
                 \end{bmatrix}:\ 
         \rho \in \mR^{r\x r},\ 
         \sigma \in \mR^{r\x m},\ 
         \tau \in \mR^{p\x r}
    \right\}.
\end{equation}
This is a linear subspace of $\mR^{(r+p)\x(r+m)}$ which inherits
the Frobenius inner product of matrices.
Let $\bPi_{r,m,p}$ denote the orthogonal projection onto
$\bGamma_{r,m,p}$ which pads the bottom-right $(p\x m)$-block of a $(r+p)\x (r+m)$-matrix with zeros:
\begin{equation}
\label{bPi}
    \bPi_{r,m,p}
    \left(
                 \begin{bmatrix}
                     \rho  &  \sigma \\
                     \tau & \varpi
                 \end{bmatrix}
    \right)
    =
        \begin{bmatrix}
            \rho  &  \sigma \\
            \tau & 0
        \end{bmatrix}.
\end{equation}
%The subscripts in $\bGamma_{r,m,p}$ and $\bPi_{r,m,p}$ will sometimes be omitted for brevity.
The dependence of the   closed-loop system matrices $\cA$, $\cB$,
$\cC$ on the controller matrices $a$, $b$, $c$ in  (\ref{F}) can be written as
\begin{equation}
\label{ABC}
    \Gamma
     :=
        \begin{bmatrix}
            \cA & \cB\\
            \cC & 0
        \end{bmatrix}
    =
    \Gamma_0
    +
    \Gamma_1
    \gamma
    \Gamma_2,
    \qquad
    \gamma
    :=
        \begin{bmatrix}
            a & b\\
            c & 0
        \end{bmatrix}.
\end{equation}
The affine map $\bGamma_{n,p_2,m_2}\ni \gamma \mapsto \Gamma \in
\bGamma_{2n,m_1,p_1}$ is specified completely  by three
matrices
\begin{equation}
\label{Gamma012}
    \Gamma_0
    :=
        \begin{bmatrix}
            \!\!0_n \, \ 0 \ \ \, 0\\
            0 \ \  A \ \  B_1\\
            \!\!\!0 \ \  C_1 \ \, 0
        \end{bmatrix},
        \qquad
    \Gamma_1
    :=
        \begin{bmatrix}
            I_n & \!\!\!0\\
            0 & \!\!\!B_2\\
            0 & \!\!\!D_{12}
        \end{bmatrix},\qquad
    \Gamma_2
    :=
        \begin{bmatrix}
            \!\!\!\!\!\!\!I_n \  \, 0 \ \  \ 0\\
            0 \ \ C_2 \  D_{21}
        \end{bmatrix},
\end{equation}
where $0_n$ denotes the $(n\x n)$-matrix of zeros.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Equations for optimal controller}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


We now obtain necessary conditions of optimality in the class
(\ref{K})  of  $n$-dimensional stabilizing controllers $K$ for the LQQG problem (\ref{LQQG}). To this end,  we compute the Frechet
derivatives of the quadro-quartic functional
 of the closed-loop system $F$
as a composite function $\gamma \mapsto \Gamma \mapsto \bQ$ of the controller matrices $a$, $b$, $c$ and equate the derivatives to zero.
The differentiation is carried out in two steps: we first consider $\cA$, $\cB$, $\cC$ to be independent variables, and then take into account their dependence on $a$, $b$, $c$.


%====================================================================
\begin{lemma}
\label{lem:dLdGamma} The Frechet derivatives of the quadro-quartic functional $\bQ$ with respect
to the closed-loop system matrices $\cA$, $\cB$, $\cC$, assembled into the matrix $\Gamma$
in (\ref{ABC}), are computed as
\begin{equation}
\label{dLdGamma}
    \d_{\Gamma} \bQ
     :=
        \begin{bmatrix}
             \d_{\cA} \bQ & \d_{\cB} \bQ \\
             \d_{\cC} \bQ & 0
        \end{bmatrix}
    =
     2
        \begin{bmatrix}
             R  & \Omega \cB\\
            \cC \Ups & 0
        \end{bmatrix}.
\end{equation}
Here,
\begin{eqnarray}
\label{Ups}
    \Ups
     & := &
     P  + \theta ( P  H  + \Phi),\\
\label{Omega}
    \Omega
     & := &
     Q  + \theta ( H  Q  + \Psi),\\
\label{R}     R  & := &
     H + \theta (H^2 + Q\Phi + \Psi P),%\\
%
%& = &
%     Q  \Ups + \Omega P  -  H  (I+\theta  H)
\end{eqnarray}
with $P$, $Q$ the Gramians from (\ref{cPQ}); the matrix $H$ is given by (\ref{H}), and $\Phi$, $\Psi$ are the controllability and observability Schattenians of $F$ satisfying the algebraic Lyapunov equations
\begin{equation}
\label{Phi_Psi}
    \cA\Phi + \Phi\cA^{\rT} +  P \cC^{\rT}\cC  P
     = 0,
     \qquad
    \cA^{\rT}\Psi
    +
    \Psi\cA +  Q \cB\cB^{\rT}  Q
     =  0.
\end{equation}
\end{lemma}
%====================================================================

\begin{proof}
By recalling (\ref{QQ}) and applying
Lemmas~\ref{lem:diffH2}, \ref{lem:diffH4} of Appendices~B, C to the closed-loop system $F$, it follows that
$$
    \d_{\Gamma}
    \bQ
    =
    \d_{\Gamma}
    (\|F\|_2^2)
    +
    \theta
    \d_{\Gamma}
    (\|F\|_4^4)/2\\
    = 2
    \begin{bmatrix}
        H & Q\cB\\
        \cC P & 0
    \end{bmatrix}
    +
    2\theta
    \begin{bmatrix}
        H^2 + Q\Phi + \Psi P  & (HQ + \Psi)\cB\\
        \cC (PH+\Phi) & 0
    \end{bmatrix}, 
$$
which, in view of the notations (\ref{Ups})--(\ref{R}), implies (\ref{dLdGamma}).
\end{proof}
%====================================================================

%By combining Lemmas~\ref{lem:PhiPsi} and \ref{lem:dLdGamma}, it
%follows that at a critical point of the Lagrange function $\Lambda$,
%the matrix $ R $ from (\ref{R}) takes the form
%\begin{multline*}
%     R
%    =
%    ( Q  + \theta( H  Q +\Psi)) P
%    +
%     Q
%    ( P  + \theta( P  H +\Phi))\\
%    -
%     H (I+\theta  H )
%    =
%     H (I+\theta  H )
%    +
%    \theta
%    (\Psi P  +  Q  \Phi).
%\end{multline*}
%==============================================================================

The Gramians $ P $, $ Q $ of the closed-loop system and related matrices (that is, $H$, $\Phi$,
$\Psi$, $\Ups$, $\Omega$, $ R $) inherit the four $(n\x n)$-block structure of the
matrix $\cA$ in (\ref{F}). The blocks are numbered as follows:
%\vskip-5mm
\begin{equation}
\label{blocks}
    \cA
    :=
    \begin{array}{cc}
    \!{}_{\leftarrow n \rightarrow}\  {}_{\leftarrow n\rightarrow} &\\
            \begin{bmatrix}
                \cA _{11} & \cA _{12}\\
                \cA _{21} & \cA _{22}
            \end{bmatrix}
    &\!\!\!\!\!
        \begin{matrix}
            \updownarrow\!{}^n\\
            \updownarrow\!{}_n
        \end{matrix}\\
        {}
    \end{array}
    =
    \begin{array}{cc}
    \!{}_{\leftarrow n \rightarrow}\  {}_{\leftarrow n\rightarrow} &\\
            \begin{bmatrix}
                \cA _{\bullet 1} & \cA _{\bullet 2}
            \end{bmatrix}
    &\!\!\!\!\!
            \updownarrow\!{}^{2n}
        \\ {}
    \end{array}
    =
    \begin{array}{cc}
    {}_{\leftarrow 2n \rightarrow}\\
            \begin{bmatrix}
                \cA _{1\bullet}\\
                \cA _{2\bullet}
            \end{bmatrix}
    &\!\!\!\!\!
        \begin{matrix}
            \updownarrow\!{}^n\\
            \updownarrow\!{}_n
        \end{matrix}\\
        {}
    \end{array}.
\end{equation}
In this notation, the $(\cdot)_{11}$ blocks are associated with the
controller state, and the $(\cdot)_{22}$ blocks pertain to the
plant state.

%====================================================================
\begin{lemma}
\label{lem:dLdgamma} The Frechet derivatives of the quadro-quartic functional $\bQ$
of the closed-loop system (\ref{F}) with respect to the
controller matrices $a$, $b$, $c$, assembled into the matrix
$\gamma$ in (\ref{ABC}), are computed as
\begin{equation}
 \label{dLdabc0}
    \d_{\gamma} \bQ
     =
        \begin{bmatrix}
             \d_a \bQ & \d_b \bQ \\
             \d_c \bQ & 0
        \end{bmatrix}
    =
    2
        \begin{bmatrix}
             R _{11}
            &  R _{12}C_2^{\rT} + \Omega_{1\bullet} \cB D_{21}^{\rT}\\
            B_2^{\rT} R _{21} +  D_{12}^{\rT}\cC \Ups_{\bullet 1} & 0
        \end{bmatrix},
\end{equation}
where the matrices $\Ups$, $\Omega$, $R$ are defined by (\ref{Ups})--(\ref{R}).
\end{lemma}
%====================================================================

\begin{proof}
 Since
$\bQ$ is a composite function of $a$, $b$, $c$ which enter this functional
 through the matrices $\cA$,
$\cB$, $\cC$ of the closed-loop system $F$, the chain rule yields
\begin{equation}
\label{dLdgamma}
    \d_{\gamma} \bQ
    =
    (
        \d_{\gamma} \Gamma
    )^{\dagger}
    (
        \d_{\Gamma} \bQ
    )
    =
    \bPi_{n,p_2,m_2}
    (
        \Gamma_1^{\rT}
        \d_{\Gamma} \bQ
        \Gamma_2^{\rT}
    ).
\end{equation}
Here, $(\cdot)^{\dagger}$ denotes the adjoint of a linear operator in the
sense of the Frobenius inner product of matrices, and
$\bPi_{n,p_2,m_2}$ is the orthogonal projection onto the subspace
$\bGamma_{n,p_2,m_2}$ defined by (\ref{bGamma})--(\ref{bPi}).
Indeed, the first variation of the affine map $\Gamma$, defined by
(\ref{ABC}), is $\delta \Gamma = \Gamma_1(\delta \gamma) \Gamma_2$.
Hence,
$
    \delta \bQ
    =
    \Tr
    (
        \d_{\Gamma} \bQ
        \delta\Gamma^{\rT}
    )
    =
    \Tr
    (
        \d_{\Gamma} \bQ
        (\Gamma_1(\delta\gamma) \Gamma_2)^{\rT}
    )
    =
    \Tr
    (
         \Gamma_1^{\rT}
        \d_{\Gamma} \bQ
        \Gamma_2^{\rT}
        \delta\gamma^{\rT}
    )
    =
    \Tr
    (
        \bPi_{n,p_2,m_2}
        (
             \Gamma_1^{\rT}
            \d_{\Gamma} \bQ
            \Gamma_2^{\rT}
        )
        \delta\gamma^{\rT}
    )
$,
which establishes (\ref{dLdgamma}). Substitution of the matrices
$\Gamma_1$ and $\Gamma_2$ from (\ref{Gamma012}) into the right-hand
side of (\ref{dLdgamma}) yields
\begin{align}
\nonumber
    \d_{\gamma} \bQ
%     =
%        \begin{bmatrix}
%             \d_a \bQ  & \d_b \bQ  \\
%             \d_c \bQ  & 0
%        \end{bmatrix}\\
    &=
    \bPi_{n,p_2,m_2}
    \left(
                 \begin{bmatrix}
                     I_n & 0 & 0\\
                     0 & B_2^{\rT} & D_{12}^{\rT}\\
                 \end{bmatrix}\!\!
                 \begin{bmatrix}
                      \d_{\cA} \bQ  & \d_{\cB} \bQ  \\
                      \d_{\cC} \bQ  & 0
                 \end{bmatrix}\!\!
                 \begin{bmatrix}
                     I_n & 0\\
                     0 & C_2^{\rT}\\
                     0 & D_{21}^{\rT}
                 \end{bmatrix}\!
    \right)\\
 \label{dLdabc}
    & =
        \begin{bmatrix}
            (\d_{\cA}\bQ )_{11} & (\d_{\cA} \bQ )_{12}C_2^{\rT} + (\d_{\cB} \bQ )_1 D_{21}^{\rT}\\
            B_2^{\rT}(\d_{\cA} \bQ )_{21} +  D_{12}^{\rT}(\d_{\cC} \bQ )_1 & 0
        \end{bmatrix}.
\end{align}
Here, in view of (\ref{dLdGamma}),
\begin{equation}
\label{dLdB1_dLdC1}
    \d_{\cA}\bQ
    = 2 R,
    \qquad
    (\d_{\cB} \bQ )_1
    =
    2\Omega_{1\bullet} \cB,
    \qquad
    (\d_{\cC} \bQ )_1
    =
    2\cC\Ups_{\bullet 1},
\end{equation}
 and the block numbering (\ref{blocks}) is used. The assertion
(\ref{dLdabc0}) of the lemma now follows from (\ref{dLdabc}) and
(\ref{dLdB1_dLdC1}).
\end{proof}
%====================================================================



 Necessary conditions for  optimality in the class of controllers (\ref{K}) for the LQQG problem (\ref{LQQG}) are now obtained by equating the blocks of the matrix $\d_{\gamma}\bQ$ in (\ref{dLdabc0})  to zero:
\begin{eqnarray}
%\label{full:cA}
%    \cA
%    & := &
%        \begin{bmatrix}
%                  a &  bC_2\\
%                  B_2 c & A
%        \end{bmatrix},\\
%\label{full:cB}
%    \cB
%    & := &
%        \begin{bmatrix}
%                  bD_{21}\\
%                  B_1
%        \end{bmatrix},\\
%\label{full:cC}
%    \cC
%    & := &
%        \begin{bmatrix}
%                  D_{12} c & C_1
%        \end{bmatrix},\\
%\nonumber\\
%\label{full:cP}
%    \cA P
%    +
%     P \cA^{\rT}
%    +
%    \cB\cB^{\rT}
%         & = & 0,\\
%\label{full:cQ}
%    \cA^{\rT} Q
%    +
%     Q  \cA
%    +
%    \cC^{\rT}\cC
%        & = & 0,\\
%\nonumber\\
%\label{full:Phit}
%    \cA\Phi + \Phi\cA^{\rT} +  P \cC^{\rT}\cC  P
%     & = &
%    0,\\
%\label{full:Psit}
%    \cA^{\rT}\Psi
%    +
%    \Psi\cA +  Q \cB\cB^{\rT}  Q
%     & = &  0,\\
%\nonumber\\
%\label{full:Ups}
%    \Ups
%      & := &
%     P  + \theta ( P  H  + \Phi),\\
%\label{full:Omega}
%    \Omega
%      & := &
%     Q  + \theta ( H  Q  + \Psi),\\
%%\nonumber\\
%\label{full:cR}
%     R  & := & H + \theta(H^2+Q  \Phi + \Psi P),\\
%\nonumber\\
\label{full:dLda}
             R _{11} = 0,\\
\label{full:dLdb}
     R _{12}C_2^{\rT} + \Omega_{1\bullet} \cB D_{21}^{\rT} = 0,\\
\label{full:dLdc}
            B_2^{\rT} R _{21} +  D_{12}^{\rT}\cC \Ups_{\bullet 1} =
            0.
\end{eqnarray}
%The transformation
%\begin{equation}
%\label{abc_transform}
%    (a,b,c)\mapsto (\sigma a\sigma^{-1},\sigma b,c\sigma^{-1}),
%\end{equation}
%with nonsingular matrices $\sigma$ of order $n$, yielding equivalent
%state-space realizations of the controller, induces the following
%symmetry group for the equations (\ref{full:cA})--(\ref{full:dLdc}):
%\begin{eqnarray}
%    (\cA,\cB,\cC)
%    & \mapsto &
%    (\Sigma\cA\Sigma^{-1}, \Sigma \cB, \cC\Sigma^{-1}),\\
%     P  &\mapsto &\Sigma P \Sigma^{\rT},\\
%     Q  &\mapsto & \Sigma^{-\rT} Q \Sigma^{-1},\\
%    \Phi &\mapsto &\Sigma\Phi\Sigma^{\rT},\\
%    \Psi &\mapsto & \Sigma^{-\rT}\Psi\Sigma^{-1},\\
%     R  &\mapsto &\Sigma^{-\rT} R \Sigma^{\rT}.
%\end{eqnarray}
%Here, $(\cdot)^{-\rT} := ((\cdot)^{-1})^{\rT}$ is the composition of
%the matrix inverse with the transpose, and
%$$
%    \Sigma
%    :=
%        \begin{bmatrix}
%            \sigma & 0\\
%            0_n & I_n
%        \end{bmatrix}.
%$$


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Observation-state and state-feedback matrices}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%==============================================================================
\begin{lemma}
\label{lem:pos_def_11} Suppose the matrix $D_{21}$ is of full
row rank, and $D_{12}$ is of full column rank. Also, let (\ref{K}) be a stabilizing controller with a minimal state-space realization. Then the top-left blocks of the matrices $ P $, $ Q $ from (\ref{cPQ}) and
$\Ups$, $\Omega$ from (\ref{Ups}), (\ref{Omega}) are all positive definite:
\begin{equation}
\label{pos_def_11}
     P _{11}\succ 0,
    \qquad
     Q _{11}\succ 0,
    \qquad
    \Ups_{11}\succ 0,
    \qquad
    \Omega_{11}\succ 0.
\end{equation}
\end{lemma}
%==============================================================================

\begin{proof}
Since $\theta\>0$,
and the matrices $ PH=PQ  P$, $HQ=Q  P Q$, associated with the Gramians $P$, $Q$, and the Schattenians $\Phi$,
$\Psi$ from (\ref{Phi_Psi}) are all positive semi-definite, then (\ref{Ups})
and (\ref{Omega}) imply that $
    \Ups\succcurlyeq  P $
    and
$
    \Omega\succcurlyeq  Q
$. Hence, the same ordering holds for the top-left blocks of these matrices:
$    \Ups_{11}\succcurlyeq  P_{11}$ and $
    \Omega_{11}\succcurlyeq  Q_{11}$.
Therefore, the last two relations in
(\ref{pos_def_11}) will follow from the first two. We will now
prove that $ P _{11}\succ 0$ under the assumptions that
$D_{21}$ is of full row rank and $(a,b)$ is controllable.
Indeed, $ P _{11}$ is the covariance matrix of the controller state:
\begin{equation}
\label{P11}
     P _{11} = \cov(\xi_t)
     =
     \frac{1}{2\pi}
     \int_{-\infty}^{+\infty}
     g(\omega)\Lambda(\omega)g(\omega)^*
     \rd \omega,
     \qquad
     g(\omega)
     :=
     (i\omega I_n-a)^{-1}b,
\end{equation}
where $
\Lambda(\omega):=h(\omega)h(\omega)^*$
 is the spectral density
associated with the observation signal $Y$ from (\ref{y}), with
$h(\omega)
    :=
    D_{21}
    +
    \begin{bmatrix}
        0 & C_2
    \end{bmatrix}
    (i\omega I_{2n}-\cA)^{-1}\cB
$.
From $\lim_{\omega\to \infty} \Lambda(\omega) = D_{21}D_{21}^{\rT}$, it follows that if $D_{21}$ is of  full row rank, then $\Lambda(\omega)\succ 0$ for all sufficiently  large $\omega$, say $|\omega|> \omega_0$. Now, if $P_{11}$ is singular, then $v^{\rT}P_{11} v = 0$ for some nonzero $v\in \mR^n$. In this case, (\ref{P11}) yields
$
    0=
    v^{\rT}
    P_{11} v
    \>
    (2\pi)^{-1}
    \int_{|\omega|>\omega_0}
    \|g(\omega)^* v\|_{\Lambda(\omega)}^2
    \rd \omega
$,
which, in view of $\Lambda(\omega)\succ 0$ over the high frequency range,   implies that $v^{\rT}g(\omega) = 0$ for all $|\omega|>\omega_0$. Hence, by considering the first $n$ terms of the Laurent series $v^{\rT} g(\omega) = \sum_{k=1}^{+\infty} v^{\rT} a^{k-1} b/(i\omega)^k$ at infinity \cite[Lemma~2.3 on pp. 16--17]{Kimura_1997}, it follows that the rank of the matrix $\left[b \mid  \ldots \mid a^{n-1}b\right]$ is less than $n$, and the pair $(a,b)$ is not controllable. Thus, the full row rank of $D_{21}$ and the controllability of $(a,b)$ indeed ensure $P_{11}\succ 0$.  By duality, a similar
reasoning shows that the observability of $(a,c)$ and the full
column rank condition on $D_{12}$ imply $ Q _{11}\succ 0$.
\end{proof}

%==============================================================================


%==============================================================================
\begin{theorem}
\label{th:bc} Suppose the matrix $D_{21}$ is of full row rank,
and $D_{12}$ is of full column rank. Then the matrices $b$ and $c$ of an optimal controller (\ref{K}) in the LQQG problem (\ref{LQQG}) with a minimal state-space
realization satisfy
\begin{eqnarray}
\label{b}
    b
    =
    -\Omega_{11}^{-1}
    (
         R _{12}C_2^{\rT}
        +
        \Omega_{12} B_1 D_{21}^{\rT}
    )
    (
        D_{21}
        D_{21}^{\rT}
    )^{-1},\\
\label{c}
    c
    =
    -
    (
        D_{12}^{\rT}
        D_{12}
    )^{-1}
    (
        B_2^{\rT}
         R _{21}
        +
        D_{12}^{\rT}C_1\Ups_{21}
    )
    \Ups_{11}^{-1},
\end{eqnarray}
where the matrices $\Ups$, $\Omega$, $R$ are defined by (\ref{Ups})--(\ref{R}).
\end{theorem}
%==============================================================================

\begin{proof}
 Substitution of
the matrices $\cB$ and $\cC$ from (\ref{F})  into (\ref{full:dLdb}) and (\ref{full:dLdc}) brings
these equations to the form
\begin{eqnarray}
\label{full:dLdb1}
     R _{12}C_2^{\rT}
    +
    (
        \Omega_{11}bD_{21}
        +
        \Omega_{12}B_1
    )
    D_{21}^{\rT}
    =
    0,\\
\label{full:dLdc1}
    B_2^{\rT} R _{21}
    +
    D_{12}^{\rT}
    (
        D_{12}c \Ups_{11}
        +
        C_1\Ups_{21}
    )
    =
    0.
\end{eqnarray}
By Lemma~\ref{lem:pos_def_11}, the matrices $\Ups_{11}$ and
$\Omega_{11}$ are nonsingular. Therefore, left multiplication of both
sides of (\ref{full:dLdb1}) by $\Omega_{11}^{-1}$ and right
multiplication by $(D_{21} D_{21}^{\rT})^{-1}$ yields (\ref{b}).
Similarly, right multiplication of both sides of (\ref{full:dLdc1})
by $\Ups_{11}^{-1}$ and left multiplication by $(D_{12}^{\rT}
D_{12})^{-1}$ yields (\ref{c}).
\end{proof}
%==============================================================================

Under the assumptions of Theorem~\ref{th:bc}, the modified set of equations for the state-space realization matrices of an optimal controller in the LQQG problem (\ref{LQQG}) is formed by the algebraic Lyapunov equations (\ref{cPQ}), (\ref{Phi_Psi}) and by the algebraic equations  (\ref{full:dLda}), (\ref{b}), (\ref{c}). In the case
 $\theta = 0$, these equations can be shown to yield the two independent Riccati equations for the standard LQG controller.

% These derivatives all vanish at a critical point of
%the Lagrange function. By equating the Frechet derivatives
%(\ref{dLdP}) and (\ref{dLdQ}) to zero, it follows that
%\begin{eqnarray}
%    \Phi
%    & = &
%     Q  + \theta  Q  P  Q  + \wt{\Phi},\\
%    \Psi
%    & = &
%     P  + \theta  P  Q  P  + \wt{\Psi},
%\end{eqnarray}
%where $\wt{\Phi}$ and $\wt{\Psi}$ are positive-semidefinite real
%symmetric  matrices satisfying the Lyapunov equations
%\begin{eqnarray}
%    \cA^{\rT}\wt{\Phi} + \wt{\Phi}\cA + 2\theta  Q \cB\cB^{\rT}  Q
%    & = & 0,\\
%    \cA\wt{\Psi} + \wt{\Psi}\cA^{\rT} + 2\theta  P \cC^{\rT}\cC  P
%    & = &
%    0.
%\end{eqnarray}
%These representations follow from the lemma below which further
%develops the Gramian anticommutation idea of Remark~2.












%M. Konstantinov, V. MehrmannCorresponding P. Petkov, ''On properties
%of Sylvester and Lyapunov operators'', Linear Algebra and its
%Applications Volume 312, Issues 1-3, 15 June 2000, Pages 35-71


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Homotopy method}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

With the matrix $\gamma$ from (\ref{ABC}), we associate a linear subspace of $\bGamma_{n,p_2,m_2}$ by
\begin{equation}
\label{mT}
    \mT(\gamma)
    =
    \left\{
        \begin{bmatrix}
            \tau a - a\tau & \tau b\\
            -c\tau & 0
        \end{bmatrix}:\,
        \tau \in \mR^{n\x n}
    \right\}.
\end{equation}
This is the tangent space  generated by the group of transformations $(a,b,c)\mapsto (\sigma a\sigma^{-1}, \sigma b, c \sigma^{-1})$ (where $\sigma \in \mR^{n\x n}$ are arbitrary nonsingular matrices),  which leave the transfer function of the controller (\ref{K}), and hence, the input-output operator of the closed-loop system (\ref{ABC}), unchanged. The matrix $\d_{\gamma}\bQ$, associated with the  controller $K$, belongs to the orthogonal complement $\mT(\gamma)^{\bot}$ of $\mT(\gamma)$ to $\bGamma_{n,p_2,m_2}$ in the sense of the Frobenius inner product. We say that the controller  delivers a strong local minimum to the quadro-quartic functional $\bQ$ in (\ref{LQQG}) if, in addition to the equality $\d_{\gamma}\bQ = 0$, it also makes the second order Frechet derivative
$
    \d_{\gamma}^2
    \bQ
    =
    \d_{\gamma}^2(\|F\|_2^2)
    +
    \theta
    \d_{\gamma}^2(\|F\|_4^4)/2
$
positive definite on the subspace $\mT(\gamma)^{\bot}$. Now, suppose there exists a smooth map $0\< \theta \mapsto \gamma_*(\theta)\in \bGamma_{n,p_2,m_2}$ such that $\gamma_*(\theta)$ is a strong local minimum  of the quadro-quartic functional $\bQ_{\theta}$ of the closed-loop system $F$ in the sense above, so that
$
    \left.\d_{\gamma}\bQ_{\theta}\right|_{\gamma = \gamma_*(\theta)} = 0
$.
By differentiating the last equality with respect to $\theta$, it follows that
\begin{equation}
\label{homo}
    \d_{\gamma}^2 \bQ_{\theta}(\gamma_*')
    +
    \d_{\gamma}(\|F\|_4^4)/2 = 0.
\end{equation}
Here, $\gamma_*'(\theta):= \d_{\theta}\gamma_*(\theta)$ and use is made of the identity $\d_{\theta} \bQ_{\theta} =\|F\|_4^4/2$ which follows from (\ref{QQ}) and, in view of the interchangeability  of the derivatives in $\theta$ and $\gamma$,  implies that $\d_{\theta}\d_{\gamma}\bQ_{\theta} = \d_{\gamma}(\|F\|_4^4)/2\in \mT(\gamma)^{\bot}$. Since the matrix $\gamma_*(\theta)$ is defined up to the orbit of the transformation group, then $\gamma_*'(\theta):= \d_{\theta}\gamma_*(\theta)$ is defined modulo the subspace $\mT(\gamma_*(\theta))$ from (\ref{mT}). Therefore, (\ref{homo}), which is a linear equation with respect to $\gamma_*'(\theta)$, can be restricted to the subspace $\mT(\gamma_*(\theta))^{\bot}$. As long as $\gamma_*(\theta)$ is a strong local minimum of $\bQ_{\theta}$, so that the self-adjoint operator $\d_{\gamma}^2\bQ$ is positive definite (and hence, invertible) on $\mT(\gamma_*(\theta))^{\bot}$, the equation  (\ref{homo}) is equivalent to
\begin{equation}
\label{homo_ODE}
    \gamma_*'(\theta)
    =
    -\bL^{-1}(\d_{\gamma}(\|F\|_4^4))/2,
\end{equation}
where $\bL$ is the restriction of $\d_{\gamma}^2\bQ$
 to the subspace $\mT(\gamma)^{\bot}$. The equation (\ref{homo_ODE}) is an ODE, with $\theta\>0$ playing the role of fictitious time. The initial value $\gamma_*(0)$ is provided by the state-space realization triple of the standard LQG controller. The computation of an LQQG controller for $\theta>0$ can be carried out by numerically integrating the homotopy ODE (\ref{homo_ODE}) initialized at $\gamma_*(0)$. The operator $\bL$ involves Frechet differentiation of solutions of algebraic Lyapunov equations with respect to their  coefficients, and the inverse $\bL^{-1}$ can be computed by using the vectorization of matrices \cite{M_1988}. The state-space formulae of the homotopy algorithm and other details of its implementation will be reported in subsequent publications.



%We also develop a homotopy method which reduces the solution
%of the set of equations to the numerical integration of the ODE
%\begin{equation}
%%\label{homo_ODE}
%    \cL_{\theta,\gamma}(\gamma\,')
%    +
%    \frac{1}{2}
%    \d_{\gamma}(\|F\|_4^4) = 0,
%    \quad
%    \gamma
%    :=
%    \left[
%        \begin{array}{cc}
%            a & b\\
%            c & 0
%        \end{array}
%    \right],
%\end{equation}
%where $\theta$ plays the role of fictitious time (with the time
%derivative $(\cdot)\,'$) and $\d_{\gamma}$ is the Frechet derivative
%with respect to the block matrix $\gamma$. Here, the second order
%Frechet derivative
%$$
%    \cL_{\theta,\gamma}
%    :=
%    \d_{\gamma}^2
%    \bQ_{\theta}(F)
%    =
%    \d_{\gamma}^2(\|F\|_2^2)
%    +
%    \frac{\theta}{2}
%    \d_{\gamma}^2(\|F\|_4^4)
%$$
%is a self-adjoint linear operator on a Hilbert space of real
%matrices with the sparsity pattern as in (\ref{homo_ODE}) and the
%inherited Frobenius inner product. With equivalent state-space
%realizations being identified, $\gamma$ delivers a strong local
%minimum to the LQQG control problem (\ref{LQQG}) if and only if
%$\cL_{\theta,\gamma}$ is essentially positive definite that can be
%monitored in the course of numerically integrating the homotopy ODE
%(\ref{homo_ODE}).




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{thebibliography}{99}
%--------------------------------------------------------------------
\bibitem{AM_1990}
B.D.O.Anderson, and J.B.Moore, {\it Optimal Control: Linear
Quadratic Methods}, Prentice Hall, Englewood Cliffs, 1990.
%%--------------------------------------------------------------------
%\bibitem{Avram_1988}
%F.Avram, On bilinear forms in Gaussian random variables and
%Toeplitz matrices, {\it Probab. Theory Relat. Fields}, vol. 79,
%1988, pp. 37--45.
%%--------------------------------------------------------------------
%\bibitem{CR_2007}
%C.D.Charalambous and  F.Rezaei, ``Stochastic uncertain systems
%subject to relative entropy constraints: induced norms and
%monotonicity properties of minimax games", {\it IEEE Trans. Autom.
%Control}, vol. 52, no. 4, pp. 647--663, Apr. 2007.
%--------------------------------------------------------------------
%\bibitem{DE_1997}
%P.Dupuis, and R.S.Ellis, {\it A Weak Convergence Approach to the
%Theory of Large Deviations}, Wiley, 1997.
%--------------------------------------------------------------------
\bibitem{DJP_2000}
P.Dupuis, M.R.James, and I.R.Petersen, Robust properties of
risk-sensitive control, {\it Math. Control Signals Systems}, vol.
13, 2000, pp. 318--332.
%--------------------------------------------------------------------
\bibitem{Frisch_1995}
U.Frisch, {\it Turbulence: The Legacy of A.N.Kolmogorov}, Cambridge
University Press, New York, 1995.
%--------------------------------------------------------------------
\bibitem{Ginovian_1994}
M.S.Ginovian, On Toeplitz type quadratic functionals of stationary
Gaussian processes, {\it Probab. Theory Relat. Fields}, vol. 100, 1994,
pp. 395--406.
%====================================================================
%\bibitem{Girsanov_1960}
%I.V.Girsanov, ``On transforming a certain class of stochastic
%processes by absolutely continuous substitution of measures'', {\it
%Theor. Probability Appl.}, vol. 5, no. 3, pp. 285--301, 1960.
%--------------------------------------------------------------------
%\bibitem{Gray_2008}
%R.M.Gray, {\it Entropy and Information Theory}, Springer-Verlag, New
%York,  2008.
%--------------------------------------------------------------------
\bibitem{GS_1958}
U.Grenander, and G.Szeg\H{o}, {\it Toeplitz Forms and Their
Applications}, University of California Press, 1958.
%--------------------------------------------------------------------
\bibitem{HJ_2007}
 R.A.Horn, and C.R.Johnson, {\it Matrix Analysis}, Cambridge
University Press, New York, 2007.
%--------------------------------------------------------------------
\bibitem{Isserlis_1918}
L.Isserlis, On a formula for the product-moment coefficient of any
order of a normal frequency distribution in any number of
variables, {\it Biometrika}, vol. 12, no. 1/2, 1918, pp. 134--139.
%--------------------------------------------------------------------
%\bibitem{Kailath_1969}
%T.Kailath, �A General Likelihood-Ratio Formula for Random Signals in
%Gaussian Noise�, {\it IEEE Trans. Information Theory}, vol. 15, no.
%3, pp. 350--361, 1969.
%--------------------------------------------------------------------
\bibitem{Kimura_1997}
H.Kimura, \textit{Chain-Scattering Approach to $\cH^{\infty}$-Control},
Birkh\"{a}user, Boston, 1997.
%--------------------------------------------------------------------
%\bibitem{Kolmogorov}
%A.N.Kolmogoroff, ``${\rm \ddot{U}}$ber die analytischen Methoden in
%der Wahrscheinlichkeitsrechnung'', {\it Ann. Math.},vol. 104, pp.
%415--458, 1931.
%--------------------------------------------------------------------
\bibitem{Magnus_1978}
J.R.Magnus,  The moments of products of quadratic forms in normal
variables, {\it Statistica Neerlandica}, vol. 32, 1978, pp. 201�210.
%==============================================================================
\bibitem{M_1988}
J.R.Magnus, \textit{Linear Structures}, Oxford University Press, New York, 1988.
%--------------------------------------------------------------------
%\bibitem{MIS_1978}
%N.N.Moiseev, Yu.P.Ivanilov, and E.M.Stolyarova, {\it Optimization
%Methods}, Nauka, Moscow, 1978.
%--------------------------------------------------------------------
\bibitem{MG_1991}
D.Mustafa, and K.Glover, {\it Minimum Entropy $\cH_{\infty}$
Control}, LNCIS, vol. 146, Springer, 1991.
%--------------------------------------------------------------------
%\bibitem{UP_2001}
% V.A.Ugrinovskii and I.R.Petersen, ``Minimax LQG
%control of stochastic partially observed uncertain systems", {\it
%SIAM J. Control Optim.}, vol. 40, no. 4, pp. 1189--1226, 2001.
%--------------------------------------------------------------------
%\bibitem{SVK_1994}
%A.V.Semyonov, I.G.Vladimirov, and A.P.Kurdyukov, ``Stochastic
%approach to $H_{\infty}$-optimization'',  {\it Proc. 33rd CDC},
%Florida, USA, December 14--16, vol.~3, pp.~2249--2250, 1994.
%%--------------------------------------------------------------------
%\bibitem{Shiryaev}
%A.N.Shiryaev, {\it Probability}, Springer, New York, 1995.
%--------------------------------------------------------------------
\bibitem{Simon_2005}
B.Simon, \textit{Trace Ideals and Their Applications}, 2nd Ed.,
American Mathematical Society, Providence, RI, 2005.
%%--------------------------------------------------------------------
%\bibitem{SIG_1998}
%R.E.Skelton, T.Iwasaki, and K.M.Grigoriadis, {\it A Unified
%Algebraic Approach to Linear Control Design}, Taylor \& Francis,
%London, 1998.
%--------------------------------------------------------------------
%\bibitem{VKS_1995}
%I.G.Vladimirov, A.P.Kurdjukov, and A.V.Semyonov, ``Anisotropy of
%signals and entropy of linear time invariant systems'',  {\it
%Doklady Akad. Nauk}, vol.~342, no.~3, pp.~583--585, 1995 (in
%Russian).
%--------------------------------------------------------------------
\bibitem{VKS_1996}
I.G.Vladimirov, A.P.Kurdyukov, and A.V.Semyonov, ``On computing the
anisotropic norm of linear discrete-time-invariant systems'',  {\it
Proceedings of the 13th IFAC World Congress}, San-Francisco,
California, USA, June 30--July 5, Vol. G, 1996, pp.~179--184.
%--------------------------------------------------------------------
%\bibitem{VDK_2006}
%I.G.Vladimirov, P.Diamond, and P.Kloeden, ``Anisotropy-based robust
%performance analysis of finite horizon linear discrete time varying
%systems''. {\it Automation and Remote Control}, no.~8,
%pp.~1265--1282, 2006.
%--------------------------------------------------------------------
%\bibitem{DVKS_2001}
%P.Diamond, I.G.Vladimirov, A.P.Kurdyukov, and A.V.Semyonov,
%Ani\-so\-tro\-py-based performance analysis of linear discrete
%time invariant control systems, {\it Int. J. Contr.}, vol.~74,
%no.~1, 2001, pp.~28--42.
%--------------------------------------------------------------------
%\bibitem{Willems_1972}
%J.C.Willems, Dissipative dynamical systems. Part I: general
%theory, {\it Arch. Rational Mech. Anal.}, vol. 45, no. 22, 1972, pp.
%321--351.
%%====================================================================
%\bibitem{Arnold}
%L.Arnold, G.C.Papanicolaou, and V.Wihstutz, ``Asymptotic analysis of
%the Lyapunov exponent and rotation number of the random oscillator
%and applications'', SIAM J. Appl. Math.  {\bf 46} (1986), no.~3,
%427--450.
%%====================================================================
%\bibitem{Bird}
%R.B.Bird, W.E.Stewart, and E.N.Lightfoot, {\it Transport Phenomena},
%2nd Ed., Wiley, New York, 2002.
%%====================================================================
%\bibitem{Fedoryuk}
%M.V.Fedoryuk, {\it Asymptotic Analysis}, Springer, Berlin, 1992.
%%====================================================================
%\bibitem{Harrison}
% P.Harrison,  {\it Quantum Wells, Wires and Dots: Theoretical and Computational
% Physics}, Wiley, New  York, 2000.
%%====================================================================
%\bibitem{Jacod}
%J.Jacod and A.N.Shiryaev, {\it Limit Theorems for Stochastic
%Processes},  , 2nd Ed., Springer, Berlin, 2003.
%%====================================================================
%\bibitem{Karatzas}
%I.Karatzas and S.E.Shreve, {\it Brownian Motion and Stochastic
%Calculus}, Springer-Verlag, New York, 1988.
%%====================================================================
%\bibitem{Khasminskii}
%R.Z.Khasminskii, ``Necessary and sufficient conditions for
%asymptotic stability of linear stochastic systems'', Teor.
%Veroyatnost. i Primenen. {\bf 12} (1967), 167--172 (in Russian)
%(English translation: Theor. Probability Appl. {\bf 12} (1967),
%144--147).
%%====================================================================
%\bibitem{Khasminskii_1969}
%R.Z.Khasminskii,  {\it Stability of Systems of Differential
%Equations with Random Parametric Excitation}, Nauka, Moscow, 1969
%(in Russian).
%
%%====================================================================
%\bibitem{Kolkunov}
%V.A.Kolkunov, ``Elimination of divergences in quantum mechanics'',
%Teor. i Matem. Fiz. {\bf 3} (1970), no.~1, 72--77 (in Russian).
%%====================================================================
%\bibitem{Landau}
%L.D.Landau and E.M.Lifshits, {\it Quantum Mechanics:
%Non-relativistic Theory}, 4th Ed.,  Nauka, Moscow, 1989  (in
%Russian).
%%%====================================================================
%%\bibitem{Laparo}
%%K.A.Laparo and G.L.Blankenship, ``Almost sure instability of a class
%%of linear stochastic systems with jump process coefficients'',
%%Lecture Notes in Mathematics 1186. Lyapunov exponents (L.Arnold and
%%V.Wihstutz, eds.),  Springer, Berlin, 1986, pp.~161--190.
%%====================================================================
%\bibitem{Oseledec}
%V.I.Oseledec, ``A multiplicative ergodic theorem. Characteristic
%Lyapunov numbers for dynamical systems'', Trans. Moscow Math. Soc.
%{\bf 19} (1969), 197--232 (Translated from Russian: Trudy Moskov.
%Mat. Obsch. {\bf 19} (1968), 179--210).
%%====================================================================
%\bibitem{Pastur}
%L.Pastur and A.Figotin, {\it Spectra of Random and Almost-Periodic
%Operators},  Springer-Verlag, Berlin, 1992.
%%====================================================================
%\bibitem{Pastur_1974}
%L.Pastur, and E.P.Fel'dman, ``Wave transmittance for a thick layer
%of a randomly inhomogeneous medium'',
% Soviet Physics JETP, {\bf 40} (1974), 241--243.
%%====================================================================
%\bibitem{Pinsky}
%M.A.Pinsky, ``Instability of the harmonic oscillator with small
%noise'', SIAM J. Appl. Math. {\bf 46} (1986), no.~3, 451--463.
%%====================================================================
%\bibitem{Razavy}
%M.Razavy, {\it Quantum Theory of Tunneling}, World Scientific,
%Singapore, 2003.
%%====================================================================
%\bibitem{Rogers}
%L.C.G.Rogers and D.Williams, {\it Diffusions, Markov Processes and
%Martingales}, Cambridge University Press, Cambridge, 2000.
%%====================================================================
%\bibitem{Roy}
%D.K.Roy, {\it Quantum Mechanical Tunnelling and Its Applications},
%World Scientific, Singapore, 1986.
%%====================================================================
%\bibitem{Sakurai}
%J.J.Sakurai, {\it Modern Quantum Mechanics}, Addison-Wesley,
%Reading, Massachusetts, 1994.
%%====================================================================
%\bibitem{Shiryaev}
%A.N.Shiryaev, {\it Probability}, Springer, New York, 1996.
%%====================================================================
%\bibitem{Vladimirov}
%V.S.Vladimirov, {\it Equations of Mathematical Physics}, Mir,
%Moscow, 1984.
%%====================================================================
%\bibitem{Wihstutz}
%V.Wihstutz,  ``Parameter dependence of the Lyapunov exponent for
%linear stochastic systems. A survey'', Lecture Notes in Mathematics
%1186. Lyapunov exponents (L.Arnold and V.Wihstutz, eds.), Springer,
%Berlin, 1986, pp.~200--215.
%%====================================================================
\end{thebibliography}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\appendix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Covariance of squared norms of Gaussian
random  vectors}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\renewcommand{\theequation}{A.\arabic{equation}}
\setcounter{equation}{0}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%==============================================================================
\begin{lemma}
\label{lem:bilinear} Let $\xi$ and $\eta$ be jointly Gaussian random
vectors with zero mean. Then the covariance of their squared
Euclidean norms is expressed in terms of the Frobenius norm of their
cross-covariance matrix by
\begin{equation}
\label{bilinear}
    \cov
    (
        |\xi|^2,
        |\eta|^2
    )
    =
    2
    \|
        \cov(\xi,\eta)
    \|^2.
\end{equation}
\end{lemma}
%==============================================================================

\begin{proof}
By applying the representation \cite{Isserlis_1918} for the mixed
 moments of Gaussian random
variables in terms of their covariances to the entries of the
vectors $\xi$ and $\eta$, it follows that
$    \bE(\xi_i^2\eta_j^2)
    =
%    \bE(
%            \xi_i
%            \xi_i
%            \eta_j
%            \eta_j
%    )
%    =
    \bE
    (
        \xi_i\xi_i
    )
    \bE
    (
        \eta_j\eta_j
    )
     +
    \bE
    (
        \xi_i\eta_j
    )
    \bE
    (
        \xi_i \eta_j
    )
     +
    \bE
    (
        \xi_i\eta_j
    )
    \bE
    (
        \xi_i\eta_j
    )
    =
    \bE(\xi_i^2)
    \bE(\eta_j^2)
    +
    2
    (
        \cov
        (
            \xi_i,
            \eta_j
        )
    )^2.
$
%$    \bE(\xi_i^2\eta_j^2)
%    =
%    \bE(
%        \underbrace
%        {
%            \xi_i
%            \xi_i
%            \eta_j
%            \eta_j
%        }_{\tiny \!1 \, \ \ 2\,\ \ 3\,\ \ 4}
%    )
%    =
%    \bE
%    (
%        \underbrace{\xi_i\xi_i}_{\tiny 1\ \ 2}
%    )
%    \bE
%    (
%        \underbrace{\eta_j\eta_j}_{\tiny 3\ \ 4}
%    )
%     +
%    \bE
%    (
%        \underbrace{\xi_i\eta_j}_{\tiny 1\ \ 3}
%    )
%    \bE
%    (
%        \underbrace{\xi_i \eta_j}_{\tiny 2\ \ 4}
%    )
%     +
%    \bE
%    (
%        \underbrace{\xi_i\eta_j}_{1\ \ 4}
%    )
%    \bE
%    (
%        \underbrace{\xi_i\eta_j}_{2\ \ 3}
%    )
%    =
%    \bE(\xi_i^2)
%    \bE(\eta_j^2)
%    +
%    2
%    (
%        \cov
%        (
%            \xi_i,
%            \eta_j
%        )
%    )^2.
%$
%\vskip-8mm\begin{multline*}
%    \bE(\xi_i^2\eta_j^2)
%    =
%    \bE(
%        \overbrace
%        {
%            \xi_i
%            \xi_i
%            \eta_j
%            \eta_j
%        }^{\tiny \!1 \, \ \ 2\,\ \ 3\,\ \ 4}
%    )\\
%    =
%    \bE
%    (
%        \underbrace{\xi_i\xi_i}_{\tiny 1\ \ 2}
%    )
%    \bE
%    (
%        \underbrace{\eta_j\eta_j}_{\tiny 3\ \ 4}
%    )
%     +
%    \bE
%    (
%        \underbrace{\xi_i\eta_j}_{\tiny 1\ \ 3}
%    )
%    \bE
%    (
%        \underbrace{\xi_i \eta_j}_{\tiny 2\ \ 4}
%    )
%     +
%    \bE
%    (
%        \underbrace{\xi_i\eta_j}_{1\ \ 4}
%    )
%    \bE
%    (
%        \underbrace{\xi_i\eta_j}_{2\ \ 3}
%    )\\
%    =
%    \bE(\xi_i^2)
%    \bE(\eta_j^2)
%    +
%    2
%    (
%        \cov
%        (
%            \xi_i,
%            \eta_j
%        )
%    )^2.
%\end{multline*}\vskip-2mm\noindent
Therefore,
\begin{equation}
\label{bilinear1}
    \bE(
        |\xi|^2
        |\eta|^2
    )
     =
    \sum_{i,j}
    \bE(\xi_i^2 \eta_j^2)
     =
    \bE(|\xi|^2) \bE(|\eta|^2)
    +
    2
    \sum_{i,j}
    (
        \cov
        (
            \xi_i,
            \eta_j
        )
    )^2,
\end{equation}
where the rightmost sum is $\|\cov(\xi,\eta)\|^2$. The relation (\ref{bilinear}) is now obtained by
substituting (\ref{bilinear1}) into
$\cov(|\xi|^2,|\eta|^2) := \bE(|\xi|^2|\eta|^2) -
\bE(|\xi|^2)\bE(|\eta|^2)$. Note that (\ref{bilinear}) can also be established by using \cite[Lemma~6.2]{Magnus_1978}.
\end{proof}
%==============================================================================


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{State space formula for Frechet derivative of $\cH_2$-norm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\renewcommand{\theequation}{B.\arabic{equation}}
\setcounter{equation}{0}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%==============================================================================
\begin{lemma}
\label{lem:diffH2}
The Frechet derivative of the squared $\cH_2$-norm $E:= \|F\|_2^2$ of the system (\ref{FABC}), with $A$ Hurwitz, is computed as
\begin{equation}
\label{dEdGamma}
    \d_{\Gamma} E
    =
    2
    \begin{bmatrix}
        H & Q B\\
        C P & 0
    \end{bmatrix},
    \qquad
    \Gamma:=
    \begin{bmatrix}
        A & B\\
        C  & 0
    \end{bmatrix}.
\end{equation}
Here, the matrix $H$ is associated by (\ref{H}) with the Gramians $P$, $Q$ from
 (\ref{PQ}).
\end{lemma}
%==============================================================================

\begin{proof}
The Frechet derivative $\d_{\Gamma} E $ inherits the block structure of the matrix $\Gamma$:
\begin{equation}
\label{dEdGamma_blocks}
    \d_{\Gamma} E
    =
    \begin{bmatrix}
        \d_A E  & \d_B E \\
        \d_C E  & 0
    \end{bmatrix}.
\end{equation}
We will now compute the blocks of this matrix.  To calculate $\d_A E $, let $B$ and $C$ be fixed. Then the first variation of $ E $ with respect to $A$ is
$
    \delta  E
    =
    \Tr(C^{\rT}C\delta P)
    =
    -\Tr((A^{\rT}Q+QA)\delta P)
    =
    -\Tr(Q(A\delta P +(\delta P)A^{\rT}))
    =
    \Tr(Q((\delta A) P +P\delta A^{\rT}))
    =
    2\Tr(H\delta A^{\rT})
$,
which implies that
\begin{equation}
\label{dEdA}
    \d_A E  = 2H.
\end{equation}
Here, use has also been made of the first variation  of the Lyapunov equation for $P$ with constant $B$ which yields
$        A\delta P + (\delta P) A^{\rT} + (\delta A)P + P\delta A^{\rT} = 0
$.
To compute $\d_B E $, we fix $A$ and $C$. Then the observability Gramian $Q$, which is a function of $A$ and $C$, is also constant, and  the first variation of $ E $ with respect to $B$ is
$
    \delta  E
    =
    \Tr(Q\delta(BB^{\rT}))
    =
    \Tr(Q((\delta B)B^{\rT}+B\delta B^{\rT}))
    =
    2\Tr(Q B\delta B^{\rT})
$,
and hence,
\begin{equation}
\label{dEdB}
    \d_B E  = 2QB.
\end{equation}
The derivative $\d_C E $ is calculated by a similar reasoning. Assuming $A$ and $B$ (and so also the controllability Gramian $P$) to be fixed, the first variation of $ E $ with respect to $C$ is
$
    \delta  E
    =
    \Tr(P\delta(C^{\rT}C))
    =
    \Tr(P((\delta C^{\rT})C+C^{\rT}\delta C))
    =
    2\Tr(C P\delta C^{\rT}),
$
which implies that
\begin{equation}
\label{dEdC}
    \d_C E  = 2C P.
\end{equation}
Substitution of (\ref{dEdA})--(\ref{dEdC}) into (\ref{dEdGamma_blocks}) yields (\ref{dEdGamma}).
\end{proof}

%==============================================================================

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\subsection{First variation of Hardy-Schatten norms}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\renewcommand{\theequation}{B.\arabic{equation}}
%\setcounter{equation}{0}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%
%
%
%
%
%%==============================================================================
%\begin{lemma}
%\label{lem:H2kvariation} The first
%variation of the $2k$th power of the $\cH_{2k}$-norm of a system $F$ is given by
%\begin{equation}
%\label{varH}
%        \delta
%        (\|F\|_{2k}^{2k})
%    =
%    \frac{k}{\pi}
%    \int_{-\infty}^{\infty}
%    \Tr
%    (
%        \wh{F}
%        G^{k-1}
%        \wh{\delta F}^*
%    )
%    \rd \omega,
%\end{equation}
%where
%\begin{equation}
%\label{G}
%    G(\omega)
%    :=
%    \wh{F}(\omega)^*\wh{F}(\omega).
%\end{equation}
%\end{lemma}
%%==============================================================================
%
%\begin{proof} For a fixed but otherwise arbitrary frequency $\omega$,
%the variation of $\Tr(G(\omega)^k)$ is
%\begin{multline}
%\label{varTrPhi}
%\delta \Tr (G^k)
%=    \sum_{j=0}^{k-1}
%    \Tr
%    (
%        G^j
%        (\delta G)
%        G^{k-1-j}
%    )
%=    k
%    \Tr
%    (
%        G^{k-1}
%        \delta G
%    )\\
%=    k
%    \Tr
%    (
%        G^{k-1}
%        (
%         \wh{\delta F}^*\wh{F}
%         +
%         \wh{F}^*\wh{\delta F}
%        )
%    ),
%\end{multline}
%where we have used the cyclic property of the trace, the equality
%$
%    \delta G
%    =
%    \wh{\delta F}^*\wh{F}
%    +
%    \wh{F}^*\wh{\delta F}
%$  which follows from (\ref{G}), and the differential
%relations of matrix calculus from \cite[p.~270]{SIG_1998}. By dividing
%both parts of (\ref{varTrPhi}) by $2\pi$ and integrating them over
%the entire frequency range, it follows that
%\begin{multline}
%\label{twoints}
%    \delta
%    (\|F\|_{2k}^{2k})
%=    \frac{k}{2\pi}
%    \left(
%        \int_{-\infty}^{+\infty}
%        \Tr
%        (
%            G^{k-1}
%             \wh{\delta F}^*\wh{F}
%        )
%        \rd \omega
%        \right.
%        \\
%        \left.
%        +
%        \int_{-\infty}^{+\infty}
%        \Tr
%        (
%            G^{k-1}
%             \wh{F}^*\wh{\delta F}
%        )
%        \rd \omega
%    \right)\\
%=    \frac{k}{\pi}
%        \int_{-\infty}^{+\infty}
%        \Tr
%        (
%            \wh{F}
%            G^{k-1}
%             \wh{\delta F}^*
%        )
%        \rd \omega
%\end{multline}
%which yields (\ref{varH}), with the finiteness of the integrals
%verified using the H\"{o}lder inequality for Schatten norms
%\cite{Simon_2005}. The first two integrals  in (\ref{twoints}) are
%equal since they reduce to each other by changing the integration
%variable
% $\omega \mapsto -\omega$. Here, we have also used the
%identities $\wh{F}(-\omega) = \overline{\wh{F}(\omega)}$ and
%$G(-\omega) = \overline{G(\omega)} = G(\omega)^{\rT}$ which
%hold for systems $F$ with real impulse response, and the invariance
%of the trace under the matrix transposition.
%\end{proof}
%%==============================================================================
%
%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\subsection{Differentiation of matrix-to-system operator}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\renewcommand{\theequation}{C.\arabic{equation}}
%\setcounter{equation}{0}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%Consider a {\it matrix-to-system} operator $\bR$ which maps a block
%matrix $\Gamma$ to an LTI system $F$ with a strictly proper rational transfer
%function, whose state-space realization matrices are the blocks of
%the matrix $\Gamma$:
%\begin{equation}
%\label{bR}
%        \begin{bmatrix}
%            A & B\\
%            C & 0
%        \end{bmatrix}
%    =:
%    \Gamma
%    \mapsto
%    F:=
%    \left[
%        \begin{array}{c|c}
%            A & B\\
%            \hline
%            C & 0
%        \end{array}
%    \right].
%\end{equation}
%Here, $A\in \mR^{n\x n}$, $B\in \mR^{n\x m}$, $C\in \mR^{p\x n}$,
%with $A$ Hurwitz. Such matrices $\Gamma$ form
%an open subset of the Hilbert space $\bGamma_{n,m,p}$ introduced in Section~\ref{sec:bGamma}. The operator $\bR$ is not
%injective, since equivalent state-space realizations are mapped to
%the same system. The nonlinear operator $\bR$ is continuous in the
%topology of pointwise convergence of transfer functions in the
%complex right half-plane. The lemma below establishes the Frechet
%differentiability of $\bR$ and computes its derivative. To this end,
%we associate with (\ref{bR}) the systems
%\begin{equation}
%\label{LR}
%    L
%    :=
%    \left[
%        \begin{array}{c|c}
%            A & I_n\\
%            \hline
%            C & 0
%        \end{array}
%    \right],
%    \quad
%    R
%    :=
%    \left[
%        \begin{array}{c|c}
%            A & B\\
%            \hline
%            I_n & 0
%        \end{array}
%    \right].
%\end{equation}
%In addition to their role in the formulation of the lemma, $L$ and
%$R$ allow  the system to be represented as
%$
%    F
%    =
%    LB
%    =
%    CR
%$.
%
%
%
%
%
%
%
%%==============================================================================
%\begin{lemma}
%\label{lem:dABC} The Frechet derivative of the asymptotically stable
%system $F$ from (\ref{bR}) with respect to the matrix $\Gamma$  is
%given by
%\begin{equation}
%\label{dABC}
%    \d_{\Gamma}
%    F(\delta\Gamma)
%    =
%        \begin{bmatrix}
%            L & I_p
%        \end{bmatrix}
%    \delta\Gamma
%        \begin{bmatrix}
%            R \\
%            I_m
%        \end{bmatrix}.
%%    =
%%    \left[
%%    \begin{array}{cc|c}
%%        A & 0 & B\\
%%        \delta A & A & \delta B\\
%%        \hline
%%        \delta C & C & 0
%%    \end{array}
%%    \right].
%\end{equation}
%\end{lemma}
%%==============================================================================
%
%\begin{proof}
%For any fixed $s\in \mC$ with $\Re s\> 0$, the transfer function
%$C(sI_n-A)^{-1}B$ of the system $F$ is differentiable with
%respect to $A$, $B$, $C$, and its first variation is computed
%as
%\begin{multline}
%\label{dABC1}
%    \delta (C(sI_n-A)^{-1}B)
%    =
%    C(sI_n-A)^{-1} (\delta A)(sI_n-A)^{-1}B\\
%     +
%    C(sI_n-A)^{-1}\delta B
%    +
%    (\delta C)(sI_n-A)^{-1}B\\
%    =
%        \begin{bmatrix}
%            C(sI_n-A)^{-1} & I_p
%        \end{bmatrix}
%        \begin{bmatrix}
%            \delta A & \delta B\\
%            \delta C & 0
%        \end{bmatrix}
%        \begin{bmatrix}
%            (sI_n-A)^{-1}B \\
%            I_m
%        \end{bmatrix}.
%\end{multline}
%Here, we use the first variation $\delta (X^{-1}) =
%-X^{-1} (\delta X) X^{-1}$ for the inverse of a nonsingular matrix
%$X$; see, for example, \cite{SIG_1998}.
% Since $C(sI_n-A)^{-1}$ and $(sI_n-A)^{-1}B$ are the transfer
%functions of the systems $L$ and $R$ from (\ref{LR}), then  (\ref{dABC1}) establishes
%(\ref{dABC}).
%\end{proof}
%%The second equality in (\ref{dABC}) is obtained from
%%the first equality in (\ref{dABC1}) by noting that
%%\begin{multline}
%%\label{dABC2}
%%C(sI_n-A)^{-1} (\delta A)(sI_n-A)^{-1}B\\
%%+    C(sI_n-A)^{-1}\delta B+
%%    (\delta C)(sI_n-A)^{-1}B
%%    =
%%    \begin{bmatrix}
%%        \delta C & C
%%    \end{bmatrix}\\
%%    \x\begin{bmatrix}
%%        (sI_n-A)^{-1} & 0\\
%%        (sI_n-A)^{-1}(\delta A)(sI_n-A)^{-1} & (sI_n-A)^{-1}
%%    \end{bmatrix}
%%    \begin{bmatrix}
%%        B \\
%%        \delta B
%%    \end{bmatrix}\\
%%    =
%%    \begin{bmatrix}
%%        \delta C & C
%%    \end{bmatrix}
%%    \left(
%%        sI_{2n}
%%        -
%%         \begin{bmatrix}
%%             A & 0\\
%%             \delta A & A
%%         \end{bmatrix}
%%    \right)^{-1}
%%    \begin{bmatrix}
%%        B \\
%%        \delta B
%%    \end{bmatrix}.
%%\end{multline}
%%Indeed, the right-hand side of (\ref{dABC2}) is the transfer function of the
%%system described by the right-hand side of (\ref{dABC}).
%%Its
%%alternative proof is based on considering the equations
%%\begin{equation}
%%\label{aux_system}
%%    \dot{\xi}_t = A\xi_t + B\omega_t,
%%    \quad
%%    \zeta_t = C\xi_t + D\omega_t,
%%\end{equation}
%%where the matrices $A$, $B$, $C$, $D$ are continuously
%%differentiable functions of an independent parameter $\tau$. Using
%%the smoothness of solutions of ODE's with respect to the parameter
%%\cite{?}, differentiation of (\ref{aux_system}) in $\tau$ yields
%%\begin{equation}
%%\label{aux_system1}
%%    \dot{\xi}_t' = A'\xi_t + A\xi_t' + B'\omega_t,
%%    \quad
%%    \zeta_t' = C'\xi_t + C\xi_t' + D'\omega_t,
%%\end{equation}
%%where $(\cdot)':= \d_{\tau}$. Therefore, $\zeta'$ is the output of
%%the following system with the $2n$-dimensional combined state
%%$(\xi,\xi')$ and input $\omega$:
%%\begin{equation}
%%\label{zetadash}
%%    \zeta'
%%    =
%%    \left[
%%    \begin{array}{cc|c}
%%        A & 0 & B\\
%%        A' & A & B'\\
%%        \hline
%%        C' & C & D'
%%    \end{array}
%%    \right]
%%    \omega.
%%\end{equation}
%%The state-space realization on the right-hand side of
%%(\ref{zetadash}) describes the Gateaux derivative of the system
%%$(A,B,C,D)$ whence (\ref{dABC}) follows.
%%\end{proof}
%%==============================================================================
%
%
%%We will now consider strictly proper systems $F$ from (\ref{bR}), whose state-space representation is specified by the
%%triplet $(A,B,C)$.
%%==============================================================================
%\begin{lemma}
%\label{lem:diffH2} For a system $F$ with the state
%space realization (\ref{bR}), where $A$ is Hurwitz,
%\begin{equation}
%\label{diffH2_Gamma}
%    \d_{\Gamma} (\|F\|_2^2)
%    =
%    2
%    \begin{bmatrix}
%        H & QB\\
%        CP & 0
%    \end{bmatrix},
%\end{equation}
%where $P$ and $Q$ are the controllability and
%observability Gramians from (\ref{PQ}), and the matrix $H$ is defined by (\ref{H}).
%\end{lemma}
%%==============================================================================
%\begin{proof}
%By combining Lemmas~\ref{lem:H2kvariation} and \ref{lem:dABC}, it follows that the first variation of
%the squared $\cH_2$-norm of the system is
%\begin{multline*}
%    \delta (\|F\|_2^2)
%    =
%    \frac{1}{\pi}
%    \int_{-\infty}^{+\infty}
%    \Tr
%    (
%        \delta \wh{F}^* \wh{F}
%    )
%    \rd
%    \omega\\
%    =
%    \frac{1}{\pi}
%    \int_{-\infty}^{+\infty}
%    \Tr
%    \left(
%        \delta \Gamma^{\rT}
%            \begin{bmatrix}
%                \wh{L}^*\\
%                I_p
%            \end{bmatrix}
%        \wh{F}
%            \begin{bmatrix}
%                \wh{R}^* & I_m
%            \end{bmatrix}
%    \right)
%    \rd
%    \omega,
%\end{multline*}
%where the systems $L$ and $R$ are defined by (\ref{LR}). Hence, the Frechet derivative of $\|F\|_2^2$ with
%respect to $\Gamma$ is
%\begin{equation}
%\label{diffH2freq}
%    \d_{\Gamma}(\|F\|_2^2)
%    =
%    \frac{1}{\pi}
%    \int_{-\infty}^{+\infty}
%    \bPi_{n,m,p}
%    \left(
%            \begin{bmatrix}
%                \wh{L}^*\\
%                I_p
%            \end{bmatrix}
%        \wh{F}
%            \begin{bmatrix}
%                \wh{R}^* & I_m
%            \end{bmatrix}
%    \right)
%    \rd
%    \omega.
%\end{equation}
%where $\bPi_{n,m,p}$ is the orthogonal projection onto the subspace $\bGamma_{n,m,p}$
%defined by (\ref{bGamma})--(\ref{bPi}). By substituting the transfer functions of $L$, $R$,
%$F$ into (\ref{diffH2freq}), it follows that
%\begin{multline}
%\label{diffH2_A}
%    \d_A(\|F\|_2^2)
%=    \frac{1}{\pi}
%    \int_{-\infty}^{+\infty}
%                \wh{L}^*
%        \wh{F}
%                \wh{R}^*
%    \rd
%    \omega\\
%=    \frac{1}{\pi}
%    \int_{-\infty}^{+\infty}
%                (i\omega I_n -A)^{-*}
%                C^{\rT}
%                C
%                (i\omega I_n -A)^{-1}\\
%                \x B
%                B^{\rT}
%                (i\omega I_n -A)^{-*}
%    \rd
%    \omega
%=    \frac{1}{\pi}
%    \int_{-\infty}^{+\infty}
%    \re^{i\omega(t-u+v)}\\
%    \x\left(
%        \int_{\mR_+^3}
%        \re^{A^{\rT}t}
%                C^{\rT}
%                C
%        \re^{Au}
%BB^{\rT}
%        \re^{A^{\rT}v}
%        \rd t \rd u\rd v
%    \right)
%    \rd
%    \omega\\
%=
%    2
%        \int_{\mR_+^2}
%        \re^{A^{\rT}t}
%                C^{\rT}
%                C
%        \re^{A(t+v)}
%                B
%                B^{\rT}
%        \re^{A^{\rT}v}
%        \rd t \rd v
%        =
%        2 QP,
%\end{multline}
%where the identity
%$\int_{-\infty}^{+\infty}\re^{i\omega t}\rd \omega = 2\pi
%\delta(t)$ is used, with $\delta$ standing for the Dirac delta-function. By
%a similar reasoning, (\ref{diffH2freq}) yields
%\begin{multline}
%\label{diffH2_B}
%    \d_B(\|F\|_2^2)
%=    \frac{1}{\pi}
%    \int_{-\infty}^{+\infty}
%                \wh{L}^*
%        \wh{F}
%    \rd
%    \omega\\
%=    \frac{1}{\pi}
%    \int_{-\infty}^{+\infty}
%                (i\omega I_n -A)^{-*}
%                C^{\rT}
%                C
%                (i\omega I_n -A)^{-1}
%    \rd
%    \omega
%    B\\
%=    \frac{1}{\pi}
%    \int_{-\infty}^{+\infty}
%    \re^{i\omega(t-u)}
%    \left(
%        \int_{\mR_+^2}
%        \re^{A^{\rT}t}
%                C^{\rT}
%                C
%        \re^{Au}
%        \rd t \rd u
%    \right)
%    \rd
%    \omega B\\
%    =
%    2QB,
%\end{multline}
%\begin{multline}
%\label{diffH2_C}
%    \d_C(\|F\|_2^2)
%=    \frac{1}{\pi}
%    \int_{-\infty}^{+\infty}
%        \wh{F}                \wh{R}^*
%    \rd
%    \omega\\
%=    \frac{1}{\pi}C
%    \int_{-\infty}^{+\infty}
%                (i\omega I_n -A)^{-1}
%                BB^{\rT}
%                (i\omega I_n -A)^{-*}
%    \rd
%    \omega
%    \\
%=    \frac{1}{\pi}C
%    \int_{-\infty}^{+\infty}
%    \re^{i\omega(u-t)}
%    \left(
%        \int_{\mR_+^2}
%        \re^{At}
%                BB^{\rT}
%        \re^{A^{\rT}u}
%        \rd t \rd u
%    \right)
%    \rd
%    \omega \\
%    =
%    2CP.
%\end{multline}
%Now, (\ref{diffH2_Gamma}) is obtained by assembling (\ref{diffH2_A})--(\ref{diffH2_C}) into (\ref{diffH2freq}).
%\end{proof}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Frechet differentiation of quartic norm in state space}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\renewcommand{\theequation}{C.\arabic{equation}}
\setcounter{equation}{0}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%==============================================================================
%The matrix $H$ in (\ref{H})  satisfies an
%algebraic Sylvester equation
%\begin{equation}
%\label{Sylvester}
%    HA^{\rT} - A^{\rT} H + QBB^{\rT} - C^{\rT}C P = 0.
%\end{equation}
%This is obtained by left-multiplying the first Lyapunov equation in
%(\ref{PQ}) by $Q$ and right-multiplying the second Lyapunov
%equation by $P$, then taking the difference and recalling (\ref{H}):
%\begin{eqnarray*}
%\nonumber
%    0
%    =
%    Q
%    (
%        A P
%        +
%        PA^{\rT}
%        +
%        BB^{\rT}
%    )
%    -
%    (
%        A^{\rT}Q
%        +
%        QA
%        +
%        C^{\rT}C
%    )
%    P\\
%    =
%    H A^{\rT}
%    -
%    A^{\rT} H
%    +
%    QBB^{\rT}
%    -
%    C^{\rT}C P.
%\end{eqnarray*}
%%The equation (\ref{Sylvester}) specifies the matrix $H$ up to the
%%commutant of the matrix $A^{\rT}$.
%%==============================================================================

%==============================================================================
\begin{lemma}
\label{lem:diffH4} The Frechet derivative of the fourth power $N:= \|F\|_4^4$ of the quartic  norm of the system (\ref{FABC}), with $A$ Hurwitz, is computed as
\begin{equation}
\label{dNdGamma}
    \d_{\Gamma} N
    =
    4
    \left[
    \begin{array}{cc}
        H^2 + Q\Phi + \Psi P & (HQ + \Psi)B\\
        C(PH + \Phi) & 0
    \end{array}
    \right].
\end{equation}
Here, the matrix $H$ is associated by (\ref{H}) with the Gramians $P$, $Q$ from
 (\ref{PQ}), and $\Phi$, $\Psi$ are the Schattenians from (\ref{Schattenian_Phi_Psi}).
\end{lemma}
%==============================================================================
\begin{proof}
We will compute the Frechet derivative of $N$ by using the representation
\begin{equation}
\label{FFF}
    N
    =
    2E_1
    =
    2E_2,
    \qquad
    E_1:= \|F_1\|_2^2,
    \qquad
    E_2:= \|F_2\|_2^2,
\end{equation}
of the $\cH_4$-norm from Lemma~\ref{lem:H4state} in terms of the squared $\cH_2$-norms  of the subsidiary systems
$
    F_1
    :=
    (A,B, B^{\rT} Q)$ and $
    F_2
    :=
    (A,PC^{\rT}, C)
$
as composite functions of the matrices $A$, $B$, $C$.
Since the controllability and observability Gramians of $F_1$ are $P$ and $\Psi$, and the controllability and observability Gramians of $F_2$ are $\Phi$ and $Q$, then application of Lemma~\ref{lem:diffH2} from Appendix~B to the systems $F_1$ and $F_2$ yields
\begin{eqnarray}
\label{diffH2F1}
    \d_{\Gamma_1}E_1
    =
    2
    \begin{bmatrix}
        \Psi P & \Psi B\\
        B^{\rT} H & 0
    \end{bmatrix},
    \qquad
    \Gamma_1
    :=
    \begin{bmatrix}
    A & B\\
    B^{\rT}Q & 0
    \end{bmatrix},\ \ \\
\label{diffH2F2}
    \hskip-5mm \d_{\Gamma_2}E_2
    =
    2
    \begin{bmatrix}
        Q \Phi & HC^{\rT}\\
        C\Phi & 0
    \end{bmatrix},
    \qquad
    \Gamma_2
    :=
    \begin{bmatrix}
    A & PC^{\rT}\\
    C & 0
    \end{bmatrix}.\ \ \,
\end{eqnarray}
Suppose the matrices $A$ and $C$ are fixed and hence, so also is $Q$. Then (\ref{diffH2F1}) implies that the first variation of  $E_1$ with respect to $B$ is
\begin{align}
\nonumber
    \delta E_1
    & =
    2
    \Tr(\Psi B \delta B^{\rT})
    +
    2\Tr(B^{\rT} H \delta (B^{\rT}Q)^{\rT})\\
\nonumber
    & =
    2
    \Tr((\Psi B +
    QH^{\rT} B)\delta B^{\rT})\\
\label{deltaH2F1}
    & =
    2
    \Tr(
        (\Psi + HQ)B
        \delta B^{\rT}
    ),
\end{align}
where the identity $Q H^{\rT}= QPQ = HQ$ has also been used. From (\ref{FFF}) and (\ref{deltaH2F1}), it follows that
\begin{equation}
\label{dH4FdB}
    \d_B N
    =
    4
        (HQ+\Psi)B.
\end{equation}
Suppose the matrices $A$ and $B$ are fixed and hence, so also is $P$. Then (\ref{diffH2F2}) implies that the first variation of  $E_2$ with respect to $C$ is
\begin{align}
\nonumber
    \delta E_2
    & =
    2
    \Tr(HC^{\rT}\delta(PC^{\rT})^{\rT})
    +
    2\Tr(C\Phi\delta C^{\rT})
    \\
\nonumber
    &=
    2
    \Tr((C\Phi +
    CH^{\rT} P)\delta C^{\rT})\\
\label{deltaH2F2}
    & =
    2
    \Tr(
        C(\Phi + PH)
        \delta C^{\rT}
    ),
\end{align}
where the identity $ H^{\rT} P= PQP = PH$ has also been used. From (\ref{FFF}) and (\ref{deltaH2F2}), it follows that
\begin{equation}
\label{dH4FdC}
    \d_C N
    =
    4
        C(PH+\Phi).
\end{equation}
Now, let $B$ and $C$ be constant. Then, in view of (\ref{diffH2F1}), the variation of $E_1$ with respect to $A$ is
\begin{equation}
\label{delta_E1}
    \delta E_1
    =
    2\Tr(\Psi P \delta A^{\rT})
    +
    2
    \Tr(B^{\rT} H \delta(B^{\rT}Q)^{\rT})
    =
    2\Tr(\Psi P \delta A^{\rT})
    +
    2
    \Tr(H^{\rT}BB^{\rT} \delta Q).
\end{equation}
The first variation of the Lyapunov equation for $Q$ in (\ref{PQ}) with $C$ constant yields
$    A^{\rT}\delta Q + (\delta Q) A + (\delta A)^{\rT}Q + Q\delta A = 0
$.
Therefore,
\begin{align}
\nonumber
    \Tr((\delta Q) B B^{\rT}H )
    & =
    -\Tr((\delta Q) (AP+PA^{\rT})H )\\
\nonumber
    &=
    -\Tr((\delta Q) APH)
    -\Tr((\delta Q) PA^{\rT}H )\\
\nonumber
    & =
    \Tr(( A^{\rT}\delta Q +(\delta A)^{\rT}Q + Q\delta A)PH)
    -\Tr((\delta Q) PA^{\rT}H )\\
\nonumber
    & =
    2\Tr(H^2 \delta A^{\rT})
    +
    \Tr(P(HA^{\rT}-A^{\rT}H)\delta Q)\\
\nonumber
    & =
    2\Tr(H^2 \delta A^{\rT})
    +
    \Tr(P(C^{\rT}C P- QBB^{\rT})\delta Q)\\
\nonumber
    & =
    2\Tr(H^2 \delta A^{\rT})
    -
    \Tr((A\Phi + \Phi A^{\rT})\delta Q)
    -
    \Tr(H^{\rT}BB^{\rT}\delta Q)\\
\nonumber
    & =
    2 \Tr((H^2+Q\Phi)\delta A^{\rT})
    -
    \Tr(H^{\rT}BB^{\rT}\delta Q)\\
\label{final}
    & =
    \Tr((H^2+Q\Phi)\delta A^{\rT}).
\end{align}
%==============================================================================
Here, we have also used the definition of the controllability Schattenian $\Phi$ in (\ref{Schattenian_Phi_Psi}), and the identity
$    HA^{\rT} - A^{\rT} H =C^{\rT}C P- QBB^{\rT}
$
which is obtained from
(\ref{PQ}) and (\ref{H}) as
$    0
    =
    Q
    (
        A P
        +
        PA^{\rT}
        +
        BB^{\rT}
    )
    -
    (
        A^{\rT}Q
        +
        QA
        +
        C^{\rT}C
    )
    P
    =
    H A^{\rT}
    -
    A^{\rT} H
    +
    QBB^{\rT}
    -
    C^{\rT}C P
$. Substitution of (\ref{final}) into (\ref{delta_E1}) yields
$\delta E_1 = 2\Tr((H^2 + Q\Phi + \Psi P)\delta A^{\rT}) $, which, in view of (\ref{FFF}), implies that
\begin{equation}
\label{dH4FdA}
    \d_A N
    =
    4(H^2 + Q\Phi + \Psi P).
\end{equation}
The representation (\ref{dNdGamma}) now follows from (\ref{dH4FdB}), (\ref{dH4FdC}) and  (\ref{dH4FdA}). \end{proof}


%The equation (\ref{Sylvester}) specifies the matrix $H$ up to the
%commutant of the matrix $A^{\rT}$.
%==============================================================================


%!!!!!!!!!!!!!!!!!!!!!!!!!!!
%
%By combining Lemmas~\ref{lem:H2kvariation} and \ref{lem:dABC}, it follows that the first variation of
%the fourth power of the $\cH_4$-norm of the system is
%\begin{multline*}
%    \delta (\|F\|_4^4)
%    =
%    \frac{2}{\pi}
%    \int_{-\infty}^{+\infty}
%    \Tr
%    (
%        \delta \wh{F}^* \wh{F} G
%    )
%    \rd
%    \omega\\
%    =
%    \frac{2}{\pi}
%    \int_{-\infty}^{+\infty}
%    \Tr
%    \left(
%        \delta \Gamma^{\rT}
%            \begin{bmatrix}
%                \wh{L}^*\\
%                I_p
%            \end{bmatrix}
%        \wh{F} G
%            \begin{bmatrix}
%                \wh{R}^* & I_m
%            \end{bmatrix}
%    \right)
%    \rd
%    \omega,
%\end{multline*}
%where the systems $L$ and $R$ are defined by (\ref{LR}) and (\ref{G}) is used. Hence, the Frechet derivative of $\|F\|_4^4$ with
%respect to $\Gamma$ is
%\begin{equation}
%\label{diffH2freq}
%    \d_{\Gamma}(\|F\|_4^4)
%    =
%    \frac{2}{\pi}
%    \int_{-\infty}^{+\infty}
%    \bPi_{n,m,p}
%    \left(
%            \begin{bmatrix}
%                \wh{L}^*\\
%                I_p
%            \end{bmatrix}
%        \wh{F} G
%            \begin{bmatrix}
%                \wh{R}^* & I_m
%            \end{bmatrix}
%    \right)
%    \rd
%    \omega.
%\end{equation}
%where $\bPi_{n,m,p}$ is the orthogonal projection onto the subspace $\bGamma_{n,m,p}$
%defined by (\ref{bGamma})--(\ref{bPi}). By substituting the transfer functions of $L$, $R$,
%$F$ into (\ref{diffH2freq}), it follows that
%\begin{multline}
%\label{diffH2_A}
%    \d_A(\|F\|_2^2)
%=    \frac{2}{\pi}
%    \int_{-\infty}^{+\infty}
%                \wh{L}^*
%        \wh{F} \wh{F}^*\wh{F}
%                \wh{R}^*
%    \rd
%    \omega\\
%=    \frac{1}{\pi}
%    \int_{-\infty}^{+\infty}
%                (i\omega I_n -A)^{-*}
%                C^{\rT}
%                C
%                (i\omega I_n -A)^{-1}\\
%                \x B
%                B^{\rT}
%                (i\omega I_n -A)^{-*}
%    \rd
%    \omega
%=    \frac{1}{\pi}
%    \int_{-\infty}^{+\infty}
%    \re^{i\omega(t-u+v)}\\
%    \x\left(
%        \int_{\mR_+^3}
%        \re^{A^{\rT}t}
%                C^{\rT}
%                C
%        \re^{Au}
%BB^{\rT}
%        \re^{A^{\rT}v}
%        \rd t \rd u\rd v
%    \right)
%    \rd
%    \omega\\
%=
%    2
%        \int_{\mR_+^2}
%        \re^{A^{\rT}t}
%                C^{\rT}
%                C
%        \re^{A(t+v)}
%                B
%                B^{\rT}
%        \re^{A^{\rT}v}
%        \rd t \rd v
%        =
%        2 QP,
%\end{multline}
%where the identity
%$\int_{-\infty}^{+\infty}\re^{i\omega t}\rd \omega = 2\pi
%\delta(t)$ is used, with $\delta$ standing for the Dirac delta-function. By
%a similar reasoning, (\ref{diffH2freq}) yields
%\begin{multline}
%\label{diffH2_B}
%    \d_B(\|F\|_2^2)
%=    \frac{1}{\pi}
%    \int_{-\infty}^{+\infty}
%                \wh{L}^*
%        \wh{F}
%    \rd
%    \omega\\
%=    \frac{1}{\pi}
%    \int_{-\infty}^{+\infty}
%                (i\omega I_n -A)^{-*}
%                C^{\rT}
%                C
%                (i\omega I_n -A)^{-1}
%    \rd
%    \omega
%    B\\
%=    \frac{1}{\pi}
%    \int_{-\infty}^{+\infty}
%    \re^{i\omega(t-u)}
%    \left(
%        \int_{\mR_+^2}
%        \re^{A^{\rT}t}
%                C^{\rT}
%                C
%        \re^{Au}
%        \rd t \rd u
%    \right)
%    \rd
%    \omega B\\
%    =
%    2QB,
%\end{multline}
%\begin{multline}
%\label{diffH2_C}
%    \d_C(\|F\|_2^2)
%=    \frac{1}{\pi}
%    \int_{-\infty}^{+\infty}
%        \wh{F}                \wh{R}^*
%    \rd
%    \omega\\
%=    \frac{1}{\pi}C
%    \int_{-\infty}^{+\infty}
%                (i\omega I_n -A)^{-1}
%                BB^{\rT}
%                (i\omega I_n -A)^{-*}
%    \rd
%    \omega
%    \\
%=    \frac{1}{\pi}C
%    \int_{-\infty}^{+\infty}
%    \re^{i\omega(u-t)}
%    \left(
%        \int_{\mR_+^2}
%        \re^{At}
%                BB^{\rT}
%        \re^{A^{\rT}u}
%        \rd t \rd u
%    \right)
%    \rd
%    \omega \\
%    =
%    2CP.
%\end{multline}
%Now, (\ref{diffH2_Gamma}) is obtained by assembling (\ref{diffH2_A})--(\ref{diffH2_C}) into (\ref{diffH2freq}).

%==============================================================================

%!!!!!!!!!!!!!!!!!!
%
%With a Hurwitz matrix $A \in \mR^{n\x n}$, we associate the {\it
%inverse Lyapunov operator} $\cL_A$  on $\mR^{n\x n}$ which maps a
%matrix $M$ to the unique solution  $N$ of the algebraic Lyapunov
%equation $A N + NA^{\rT} + M = 0$:
%\begin{equation}
%\label{ILO}
%    \cL_A(M)
%    :=
%    \int_{0}^{+\infty}
%    \re^{At} M \re^{A^{\rT}t}
%    \rd t.
%\end{equation}
%The adjoint of $\cL_A$ in $\mR^{n\x n}$ (as a Hilbert space with the
%Frobenius inner product) is $\cL_A^{\dagger} = \cL_{A^{\rT}}$. Since
%$\cL_A$ commutes with the matrix transpose, that is, $\cL_A(M^{\rT})
%= (\cL_A(M))^{\rT}$, then it also does with the {\it symmetrizer}
%$\cS$ defined by
%\begin{equation}
%\label{cS}
%    \cS(M)
%    :=
%    (M+M^{\rT})/2.
%\end{equation}
%The operator $\cS$ is the orthogonal projection onto the subspace of symmetric matrices.
%The first variations of the controllability and observability
%Gramians $P$ and $Q$ of the system $(A,B,C)$ can be written using
%(\ref{ILO}) and (\ref{cS}) as
%\begin{eqnarray}
%\label{dP}
%    \delta P
%    & = &
%    2(\cL_A\circ \cS)
%    \left(
%        \left[
%            \begin{array}{cc}
%                \delta A & \delta B
%            \end{array}
%        \right]
%        \left[
%            \begin{array}{c}
%                P \\
%                B^{\rT}
%            \end{array}
%        \right]
%    \right),\\
%%    =
%%    \int_{0}^{+\infty}
%%    \re^{At}
%%    \left(
%%        \left[
%%            \begin{array}{cc}
%%                \delta A & \delta B
%%            \end{array}
%%        \right]
%%        \left[
%%            \begin{array}{c}
%%                P \\
%%                B^{\rT}
%%            \end{array}
%%        \right]
%%        +
%%        \left[
%%            \begin{array}{cc}
%%                P & B
%%            \end{array}
%%        \right]
%%        \left[
%%            \begin{array}{c}
%%                \delta A^{\rT} \\
%%                \delta B^{\rT}
%%            \end{array}
%%        \right]
%%    \right)
%%    \re^{A^{\rT}t}
%%    \rd t,\\
%\label{dQ}
%    \delta Q
%    & = &
%    2(\cL_{A^{\rT}}\circ \cS)
%    \left(
%        \left[
%            \begin{array}{cc}
%                Q & C^{\rT}
%            \end{array}
%        \right]
%        \left[
%            \begin{array}{c}
%                \delta A \\
%                \delta C
%            \end{array}
%        \right]
%    \right),
%%    \int_{0}^{+\infty}
%%    \re^{A^{\rT}t}
%%    \left(
%%        \left[
%%            \begin{array}{cc}
%%                \delta A^{\rT} & \delta C^{\rT}
%%            \end{array}
%%        \right]
%%        \left[
%%            \begin{array}{c}
%%                Q \\
%%                C
%%            \end{array}
%%        \right]
%%        +
%%        \left[
%%            \begin{array}{cc}
%%                Q & C^{\rT}
%%            \end{array}
%%        \right]
%%        \left[
%%            \begin{array}{c}
%%                \delta A \\
%%                \delta C
%%            \end{array}
%%        \right]
%%    \right)
%%    \re^{At}
%%    \rd t.
%\end{eqnarray}
%where $\circ$ is the composition of maps. Indeed,  (\ref{dP}) is
%obtained by considering the first variation of the algebraic
%Lyapunov equation for $P$ in (\ref{PQ}) which gives
%\begin{eqnarray*}
%    0
%    & = &
%    A\delta P + (\delta A)P
%    +
%    (\delta P) A^{\rT}
%     + P \delta A^{\rT}
%    +
%    (\delta B) B^{\rT} + B\delta B^{\rT}\\
%    & = &
%    A\delta P    +    (\delta P) A^{\rT}
%    +
%    \left[
%        \begin{array}{cc}
%            \delta A & \delta B
%        \end{array}
%    \right]
%    \left[
%        \begin{array}{cc}
%            P \\
%            B^{\rT}
%        \end{array}
%    \right]
%    +
%    \left[
%        \begin{array}{cc}
%            P & B
%        \end{array}
%    \right]
%    \left[
%        \begin{array}{c}
%            \delta A^{\rT} \\
%            \delta B^{\rT}
%        \end{array}
%    \right]\\
%    & = &
%    A\delta P
%    +
%    (\delta P) A^{\rT}
%    +
%    2\cS
%    \left(
%        \left[
%            \begin{array}{cc}
%                \delta A & \delta B
%            \end{array}
%        \right]
%        \left[
%            \begin{array}{c}
%                P \\
%                B^{\rT}
%            \end{array}
%        \right]
%    \right).
%\end{eqnarray*}
%This is an algebraic Lyapunov equation for $\delta P$ with the same
%matrix $A$, whence (\ref{dP}) follows. The dual relation (\ref{dQ})
%is established by a similar reasoning.

%Considering the usual block matrices,
%\begin{eqnarray*}
%    \left[
%    \begin{array}{ccc}
%        A & 0_{n\x n} & B\\
%        \delta A & A & \delta B\\
%        \delta C & C & \delta D
%    \end{array}
%    \right]
%    & = &
%    \left[
%    \begin{array}{ccc}
%        A & 0_{n\x n} & B\\
%        0_{n\x n} & A & 0_{n\x m}\\
%        0_{p\x n} & C & 0_{p\x m}
%    \end{array}
%    \right]\\
%    & + &
%    \left[
%    \begin{array}{cc}
%        0_{n\x n} & 0_{n\x p}\\
%        I_n & 0_{n\x p} \\
%        0_{p\x n} & I_p
%    \end{array}
%    \right]
%    \left[
%    \begin{array}{cc}
%        \delta A  & \delta B\\
%        \delta C & \delta D
%    \end{array}
%    \right]
%    \left[
%    \begin{array}{ccc}
%        I_n & 0_{n\x n} & 0_{n\x m}\\
%        0_{m\x n} & 0_{m\x n} & I_m
%    \end{array}
%    \right].
%\end{eqnarray*}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section{Standard LQG control problem revisited}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\setcounter{equation}{0}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%%$$
%%     R  :=  Q  P
%%$$
%
%$$
%     Q _{1\bullet}  P _{\bullet 1}
%    =
%     Q _{11}  P _{11}
%    +
%     Q _{12}  P _{21}= 0
%$$
%
%$$
%     Q _{1\bullet}( P _{\bullet 2} C_2^{\rT} + \cB D_{21}^{\rT}) = 0
%$$
%
%$$
%    (B_2^{\rT} Q _{2\bullet}  + D_{12}^{\rT}\cC ) P _{\bullet 1} = 0
%$$
%
%
%$$
%    b
%    +
%    (
%         P _{12} C_2^{\rT}
%        +
%        \Sigma
%        (
%             P _{22}C_2^{\rT}
%            +
%            B_1 D_{21}^{\rT}
%        )
%    )
%    (
%        D_{21} D_{21}^{\rT}
%    )^{-1}
%    =0
%$$
%
%
%$$
%    c
%    +
%    (
%        D_{12}^{\rT}D_{12}
%    )^{-1}
%    (
%        B_2^{\rT} Q _{21}
%        +
%        (
%            B_2^{\rT}  Q _{22}
%            +
%            D_{12}^{\rT} C_1
%        )
%        \Theta
%    )
%    =0
%$$
%
%\fbox{ $
%    \Sigma
%    :=
%     Q _{11}^{-1}
%     Q _{12},
%    \quad
%    \Theta
%    :=
%     P _{21} P _{11}^{-1}
%$}
%
%\fbox{$
%    I + \Sigma \Theta = 0
%$}
%
% Schur complements of the top-left blocks of the Gramians $ P $
%and $ Q $:
%\begin{eqnarray}
%     P _{2|1}
%& := &
%     P _{22}- P _{21} P _{11}^{-1} P _{12}
%    =
%     P _{22}-\Theta P _{12},\\
%\qquad
%     Q _{2|1}
%& := &
%     Q _{22}- Q _{21} Q _{11}^{-1} Q _{12}
%    =
%     Q _{22}- Q _{21}\Sigma
%\end{eqnarray}
%
%\fbox{$
%    b
%    +
%        \Sigma
%        (
%             P _{2|1}C_2^{\rT}
%            +
%            B_1 D_{21}^{\rT}
%        )
%    (
%        D_{21} D_{21}^{\rT}
%    )^{-1}
%    =0
%$}
%
%
%\fbox{$
%    c
%    +
%    (
%        D_{12}^{\rT}D_{12}
%    )^{-1}
%        (
%            B_2^{\rT}  Q _{2|1}
%            +
%            D_{12}^{\rT} C_1
%        )
%        \Theta
%    =0
%$}
%
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%
% We denote by $\fY_t$ the
%$\sigma$-algebra of random events generated by the past history of
%the observation signal $Y$ up until time $t$. The filtering estimate
%of the state $x_t$ of the system at time $t$ is
%$$
%    \wh{x}_t
%    :=
%    \bE(x_t\mid \fY_t).
%$$
%The process $\wh{X}$ can be represented as the output of a state
%estimator $E$ with inputs $Y$ and $U$ governed by the SDE
%$$
%    \rd \wh{x}_t
%    =
%    (A\wh{x}_t + B_2 u_t)\rd t
%    +
%    \wh{b}(\rd y_t - C_2 \wh{x}_t \rd t).
%$$
%where
%$$
%    \wh{b}
%    :=
%    (\Pi C_2^{\rT} + B_1D_{21}^{\rT})
%    (D_{21}D_{21}^{\rT})^{-1}
%$$
%and $\Pi$ is a positive semi-definite matrix which is a unique
%stabilizing solution of the algebraic Riccati equation
%$$
%    A\Pi+\Pi A^{\rT} + B_1B_1^{\rT}
%    -
%    (\Pi C_2^{\rT} + B_1D_{21}^{\rT})
%    (D_{21}D_{21}^{\rT})^{-1}
%    (\Pi C_2^{\rT} + B_1D_{21}^{\rT})^{\rT} = 0.
%$$
%in the sense that $A-\wh{b}C_2$ is Hurwitz.
%
%$$
%    e_t
%    :=
%    x_t - \wh{x}_t
%$$
%
%$$
%    \rd e_t
%    =
%    (A-bC_2)e_t\rd t
%    +
%    (B_1-bD_{21})\rd w_t
%$$
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section{Passage to Schur complements}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%\begin{eqnarray*}
%    P
%    & := &
%    \left[
%        \begin{array}{cc}
%            P_{11} & P_{12}\\
%            P_{21} & P_{22}
%        \end{array}
%    \right]\\
%    & = &
%    \left[
%        \begin{array}{cc}
%            I & 0\\
%            \Theta & I
%        \end{array}
%    \right]
%    \left[
%        \begin{array}{cc}
%            P_1 & 0\\
%            0 & P_2
%        \end{array}
%    \right]
%    \left[
%        \begin{array}{cc}
%            I & \Theta^{\rT}\\
%            0 & I
%        \end{array}
%    \right]
%     =
%    \left[
%        \begin{array}{cc}
%            P_1 & P_1\Theta^{\rT}\\
%            \Theta P_1 & \Theta P_1 \Theta ^{\rT}+P_2
%        \end{array}
%    \right]
%\end{eqnarray*}
%\begin{eqnarray*}
%    R
%    & := &
%    \left[
%        \begin{array}{cc}
%            R_{11} & R_{12}\\
%            R_{21} & R_{22}
%        \end{array}
%    \right]\\
%    & = &
%    \left[
%        \begin{array}{cc}
%            I & 0\\
%            \Delta & I
%        \end{array}
%    \right]
%    \left[
%        \begin{array}{cc}
%            R_1 & 0\\
%            0 & R_2
%        \end{array}
%    \right]
%    \left[
%        \begin{array}{cc}
%            I & \Delta^{\rT}\\
%            0 & I
%        \end{array}
%    \right]
%     =
%    \left[
%        \begin{array}{cc}
%            R_1 & R_1\Delta^{\rT}\\
%            \Delta R_1 & \Delta R_1 \Delta^{\rT}+R_2
%        \end{array}
%    \right]
%\end{eqnarray*}
%
%
%\begin{eqnarray*}
%    \alpha
%    & := &
%    \left[
%        \begin{array}{cc}
%            \alpha_{11} & \alpha_{12}\\
%            \alpha_{21} & \alpha_{22}
%        \end{array}
%    \right]\\
%    & = &
%    \cA
%    \left[
%        \begin{array}{cc}
%            I & 0\\
%            \Theta & I
%        \end{array}
%    \right]
%    =
%    \left[
%        \begin{array}{cc}
%            a+bC_2\Theta & bC_2\\
%            B_2c + A\Theta & A
%        \end{array}
%    \right]
%\end{eqnarray*}
%\begin{eqnarray*}
%    \beta
%    & := &
%    \left[
%        \begin{array}{c}
%            \beta_1 \\
%            \beta_2
%        \end{array}
%    \right]\\
%    & = &
%    \left[
%        \begin{array}{cc}
%            I & 0\\
%            -\Theta & I
%        \end{array}
%    \right]
%    \cB
%    =
%    \left[
%        \begin{array}{cc}
%            bD_{21}\\
%            B_1-\Theta bD_{21}
%        \end{array}
%    \right]
%\end{eqnarray*}
%
%
%$$
%    \Theta := P_{21}P_{11}^{-1},
%    \qquad
%    P_1:= P_{11},
%    \qquad
%    P_2
%    :=
%    P_{22}-P_{21}P_{11}^{-1}P_{12}
%    =
%    P_{22}-\Theta P_1\Theta^{\rT}
%$$
%
%\begin{eqnarray}
%    \dot{P}_1
%    & = &
%    \alpha_{11} P_1
%    +
%    P_1
%    \alpha_{11}^{\rT}
%    +
%    bD_{21}D_{21}^{\rT}b^{\rT},\\
%    \dot{P}_{21}
%    & = &
%    \alpha_{21}P_1
%    +
%    \Theta P_1 \alpha_{11}^{\rT}
%    +
%    (P_2 C_2^{\rT}
%    +
%    B_1D_{21}^{\rT})
%    b^{\rT}
%\end{eqnarray}
%
%%$$
%%    \dot{P}_{21}
%%    =
%%    (B_2c + A\Theta)P_1
%%    +
%%    \Theta P_1 a^{\rT}
%%    +
%%    (P_2 + \Theta P_1 \Theta^{\rT})C_2^{\rT}b^{\rT}
%%    +
%%    B_1D_{21}^{\rT}b^{\rT}
%%$$
%
%\begin{eqnarray*}
%    \dot{\Theta}
%    P_1
%    & = &
%    \dot{P}_{21}
%    -
%    \Theta \dot{P}_1\\
%    & = &
%    (B_2c+A\Theta -\Theta (a+bC_2\Theta))P_1
%    +
%    (
%    P_2C_2^{\rT}
%    +
%    (B_1-\Theta b D_{21})D_{21}^{\rT}
%    )
%    b^{\rT}
%\end{eqnarray*}
%
%$$
%    \dot{P}_{22}
%    =
%    (B_2c+A\Theta)P_1\Theta^{\rT}
%    +
%    \Theta P_1(B_2c+A\Theta)^{\rT}
%    +
%    AP_2
%    +
%    P_2 A^{\rT}
%    +
%    B_1B_1^{\rT}
%$$
%
%\begin{eqnarray*}
%    \dot{P}_2
%    & = &
%    \dot{P}_{22}
%    -\dot{\Theta}P_1 \Theta^{\rT}
%    -\Theta P_1 \dot{\Theta}^{\rT}
%    -\Theta \dot{P}_1 \Theta^{\rT}\\
%    & = &
%    (B_2c+A\Theta)P_1\Theta^{\rT}
%    +
%    \Theta P_1(B_2c+A\Theta)^{\rT}
%    +
%    AP_2
%    +
%    P_2 A^{\rT}
%    +
%    B_1B_1^{\rT}\\
%    & - &
%    (B_2c+A\Theta)P_1\Theta^{\rT}
%    +
%    \Theta (a+bC_2\Theta)P_1\Theta^{\rT}
%    -
%    (
%    P_2C_2^{\rT}
%    +
%    (B_1-\Theta b D_{21})D_{21}^{\rT}
%    )
%    b^{\rT}\Theta^{\rT}\\
%& - &
%    \Theta P_1(B_2c+A\Theta)^{\rT}
%    +
%    \Theta P_1(a+bC_2\Theta)^{\rT}\Theta^{\rT}
%    -
%    \Theta b
%    (
%    P_2C_2^{\rT}
%    +
%    (B_1-\Theta b D_{21})D_{21}^{\rT}
%    )^{\rT}\\
%& - &
%    \Theta(    (a+bC_2\Theta)P_1
%    +
%    P_1
%    (a+bC_2\Theta)^{\rT}
%    +
%    b
%    D_{21}D_{21}^{\rT}b^{\rT}
%)\Theta^{\rT}\\
%& = &
%    AP_2
%    +
%    P_2 A^{\rT}
%    +
%    B_1B_1^{\rT}
%    -(P_2C_2^{\rT} + B_1D_{21}^{\rT})b^{\rT}\Theta^{\rT}
%    -\Theta b(P_2C_2^{\rT} + B_1D_{21}^{\rT})^{\rT}
%    +
%    \Theta bD_{21}D_{21}^{\rT}b^{\rT}\Theta^{\rT}\\
%& = &
%    AP_2
%    +
%    P_2 A^{\rT}
%    +
%    B_1B_1^{\rT}
%    -
%    \Theta
%        \wh{b}
%        D_{21}D_{21}^{\rT}
%        \wh{b}^{\rT}
%    \Theta^{\rT}
%     +
%    \Theta
%        (b-\wh{b})
%        D_{21}D_{21}^{\rT}
%        (b-\wh{b})^{\rT}
%    \Theta^{\rT}
%\end{eqnarray*}
%where
%$$
%    \wh{b}
%    :=
%    \Theta^{-1}
%    (P_2C_2^{\rT} + B_1D_{21}^{\rT})
%    (D_{21}D_{21}^{\rT})^{-1}
%$$
%is the gain matrix in the Kalman filter
%
%
%
%




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section{Finite-horizon version of LQQG problem}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%We will now consider a variant of the LQQG control on a bounded time
%interval $[0,T]$ for a fixed but otherwise arbitrary finite time
%horizon $T>0$. This setting is formulated as
%\begin{equation}
%\label{QQT}
%    \bQ_{\theta,T}
%    :=
%    \bE\cE_T + \frac{\theta}{4}\var(\cE_T)
%    \longrightarrow
%    \min.
%\end{equation}
%Here, the minimization  is carried out over {\it nonrandom}
%functions of time $a_t$, $b_t$ and $c_t$, which take values in
%$\mR^{n\x n}$, $\mR^{n\x p_2}$ and $\mR^{m_2\x n}$, respectively,
%and specify the dynamics of the controller by
%\begin{equation}
%\label{LTV_controller}
%    \rd \xi_t
%     =
%    a_t\xi_t\rd t + b_t\rd y_t,
%    \qquad
%    u_t
%     =
%    c_t\xi_t.
%\end{equation}
%These equations are similar to  (\ref{y_to_xi_yxi_to_u}) except that
%the controller, considered in the finite-horizon setting
%(\ref{QQT}), is time-varying. The state-space realization matrices
%of the closed-loop system (\ref{F}) also become time-dependent:
%\begin{equation}
%\label{ABCt}
%    \left[
%        \begin{array}{c|c}
%            \cA_t & \cB_t\\
%            \hline
%            \cC_t & 0
%        \end{array}
%    \right]
%    =
%    \left[
%        \begin{array}{cc|c}
%                  a_t &  b_tC_2 & b_tD_{21}\\
%                  B_2 c_t & A & B_1\\
%                  \hline
%                  D_{12} c_t & C_1 & 0
%        \end{array}
%    \right].
%\end{equation}
%The combined state
%$$
%    \cX_t
%    :=
%    \left[
%        \begin{array}{c}
%          \xi_t \\
%          x_t
%        \end{array}
%    \right]
%$$
%of the closed-loop system is a Gaussian Markov process driven by the
%linear Ito SDE with time-varying coefficients and the output $Z$ is
%a time-varying linear memoryless function of the state $\cX$:
%\begin{eqnarray}
%\label{cX_SDE}
%    \rd\cX_t
%    & = &
%    \cA_t\cX_t\rd t
%    +
%    \cB_t \rd W_t,\\
%\label{z_cX}
%    \qquad
%    z_t
%    & = &
%    \cC_t \cX_t.
%\end{eqnarray}
%The covariance function of the state of the closed-loop system
%\begin{equation}
%\label{cPst}
%     P _{s,t}
%    :=
%    \cov(\cX_s, \cX_t)
%    =
%     P _{t,s}^{\rT},
%\end{equation}
%satisfies the PDE's
%\begin{equation}
%\label{cP_PDE}
%    \d_t P _{t,s}
%    =
%    \cA_t P _{t,s},
%    \quad
%    \d_t P _{s,t}
%    =
%     P _{s,t}\cA_t^{\rT},
%    \qquad
%    s< t.
%\end{equation}
%Here, the second PDE is obtained by transposing the first one in
%view of the rightmost equality from (\ref{cPst}). The state
%covariance matrix
%\begin{equation}
%\label{PP}
%     P _t:=  P _{t,t},
%\end{equation}
%which specifies the common boundary condition for the PDE's
%(\ref{cP_PDE}), is governed by the differential Lyapunov equation
%\begin{equation}
%\label{dotP}
%    \dot{ P }_t
%    =
%    \cA_t P _t
%    +
%     P _t\cA_t^{\rT}
%    +
%    \cB_t\cB_t^{\rT}.
%\end{equation}
%
%%==============================================================================
%\begin{lemma}
%\label{lem:QQrate} The quadro-quartic functional $\bQ_{\theta,T}$,
%defined by (\ref{QQT}), is representable as
%\begin{equation}
%\label{QQrate}
%    \bQ_{\theta,T}
%    =
%    \int_{0}^{T}
%    \Tr
%    (
%        \cC_t
%        ( P _t + \theta  R _t)
%        \cC_t^{\rT}
%    )
%    \rd t,
%\end{equation}
%where $ R _t$ is a real positive semi-definite symmetric
%matrix-valued function satisfying the ODE
%\begin{equation}
%\label{dotR}
%    \dot{ R }_t
%    =
%    \cA_t R _t
%    +
%     R _t\cA_t^{\rT}
%    +
%     P _t\cC_t^{\rT}\cC_t P _t,
%    \qquad
%     R _0=0.
%\end{equation}
%\end{lemma}
%%==============================================================================
%
%\textit{Proof of Lemma~\ref{lem:QQrate}:} Recall that
%$(\var(\zeta))^{^\bullet} = 2\cov(\zeta, \dot{\zeta})$ holds for any
%mean square differentiable random process $\zeta$; see, for example,
%\cite{?}. Applying this differential relation to the output energy
%process $\cE_t$ and using $\dot{\cE}_t = |z_t|^2$ (which follows
%from (\ref{ET})) and the bilinearity of the covariance operator and
%Lemma~\ref{lem:bilinear} of Appendix~A, it follows that
%%the output
%%energy variance production rate can be computed as
%\begin{eqnarray}
%\nonumber
%    (\var(\cE_t))^{^{\bullet}}
%& =&
%    2\cov(\cE_t, \dot{\cE}_t)\\
%\nonumber
%    & = &
%    2
%    \int_{0}^{t}
%    \cov(|z_s|^2,|z_t|^2)
%    \rd s\\
%\nonumber
%& = &
%    4
%    \int_{0}^{t}
%    \|
%        \cov(z_s,z_t)
%    \|^2
%    \rd s\\
%\label{dotvarE}
%& = &
%    4
%    \int_{0}^{t}
%    \Tr
%    (
%        \cC_t
%         P _{t,s}
%        \cC_s^{\rT}
%        \cC_s
%         P _{s,t}
%        \cC_t
%    )
%    \rd s
%     =
%    4
%    \Tr
%    (
%        \cC_t R _t\cC_t^{\rT}
%    ).
%\end{eqnarray}
%Here, we have also used the relation $\cov(z_s,z_t) = \cC_s
% P _{s,t} \cC_t^{\rT}$,  which follows from (\ref{z_cX}), where
%$ P _{s,t}$ is the state covariance function of the closed-loop
%system, defined by (\ref{cPst}), and
%\begin{equation}
%\label{cRt}
%     R _t
%    :=
%    \int_{0}^{t}
%         P _{t,s}
%        \cC_s^{\rT}
%        \cC_s
%         P _{s,t}
%        \rd s.
%\end{equation}
%The properties $ R _0=0$ and $ R _t =  R _t^{\rT}\succcurlyeq 0$
%follow immediately from the definition. By differentiating
%(\ref{cRt}) with respect to $t$ and using the PDE's (\ref{cP_PDE})
%together with (\ref{PP}), it follows that
%\begin{eqnarray}
%\nonumber
%    \dot{R}_t
%    & = &
%    \int_{0}^{t}
%        (\d_t  P _{t,s})
%        \cC_s^{\rT}
%        \cC_s
%         P _{s,t}
%        \rd s\\
%    \nonumber
%        & + &
%    \int_{0}^{t}
%         P _{t,s}
%        \cC_s^{\rT}
%        \cC_s
%        \d_t  P _{s,t}
%        \rd s\\
%        & + &
%         P _{t,t}
%        \cC_t^{\rT}
%        \cC_t
%         P _{t,t}
%         =
%    \cA_t R _t
%    +
%     R _t\cA_t^{\rT}
%    +
%     P _t
%    \cC_t^{\rT}
%    \cC_t
%     P _t,
%\end{eqnarray}
%which replicates the ODE in (\ref{dotR}). Substitution of
%$\bE(|z_t|^2) = \Tr(\cC_t P _t\cC_t^{\rT})$ and (\ref{dotvarE}) into
%the quadro-quartic functional from (\ref{QQT}) yields
%\begin{eqnarray*}
%    \bQ_{\theta,T}
%    & = &
%    \int_{0}^{T}
%    \left(
%        \bE(|z_t|^2)
%        +
%        \frac{\theta}{4}
%        (\var(\cE_t))^{^{\bullet}}
%    \right)
%    \rd t\\
%    & = &
%    \int_{0}^{T}
%    \Tr
%    (
%        \cC_t
%        (
%             P _t + \theta R _t
%        )
%        \cC_t^{\rT}
%    )
%    \rd t,
%\end{eqnarray*}
%which establishes (\ref{QQrate}).
%\endproof
%%==============================================================================
%
%As the proof of Lemma~\ref{lem:QQrate} shows, the quantity     $4
%    \Tr
%    (
%        \cC_t R _t\cC_t^{\rT}
%    )$ on the right-hand side of (\ref{dotvarE}) describes the output
%    energy variance production rate of the system (\ref{ABCt}) at time
%    $t$. Furthermore,  (\ref{dotP}) and (\ref{dotR}) imply that $ P _t + \theta R _t$ is the
%controllability Gramian of a subsidiary system
%$$
%        \left[
%            \begin{array}{c|cc}
%            \cA_t & \cB_t & \sqrt{\theta}  P _t\cC_t^{\rT}\\
%            \hline
%            \cC_t & 0 & 0
%            \end{array}
%        \right],
%$$
%%Similarly, the block-diagonal matrix $\diag( P _t,  R _t)$ is the
%%controllability Gramian of yet another system
%%$$
%%        \left[
%%            \begin{array}{cc|cc}
%%            \cA_t & 0 & \cB_t & 0\\
%%            0 & \cA_t & 0 &  P _t\cC_t^{\rT}\\
%%            \hline
%%            \cC_t & \sqrt{\theta} \cC_t & 0 & 0
%%            \end{array}
%%        \right]
%%$$
%%whose LQG cost  coincides with $\bQ_{\theta,T}$.
%and the expected value of the output energy of this system over the
%time interval $[0,T]$ coincides with the quadro-quartic functional
%$\bQ_{\theta,T}$ associated with the closed-loop system
%(\ref{ABCt}). The minimal value of this functional in the LQQG
%problem (\ref{QQT}) is
%$$
%    \inf \bQ_{\theta,T}
%    =
%    J_T( P _0, 0)
%$$
%where
%$$
%     R _0
%    =
%    \left[
%    \begin{array}{cc}
%        0 & 0\\
%        0 & \Pi_0
%    \end{array}
%    \right].
%$$
%Here, $\Pi_0$ is the the prior covariance matrix of the plant state
%and the controller (\ref{LTV_controller}) is assumed to be
%zero-initialized $\xi_0= 0$.
%
%
%
%\begin{equation}
%\label{JT}
%    J_T(P, R)
%    :=
%    \inf_{a_{\bullet},b_{\bullet},c_{\bullet}}
%    \int_{0}^{T}
%    \Tr
%    (
%        \cC_t
%        ( P _t + \theta  R _t)
%        \cC_t^{\rT}
%    )
%    \rd t,
%\end{equation} where the minimization is carried out over the
%functions $a_t$, $b_t$ and $c_t$ on the time interval $0\< t\< T$,
%which specify the right-hand sides of the ODE's
%\begin{eqnarray}
%\label{dotP1}
%    \dot{ P }_t
%    & = &
%    \cA_t P _t
%    +
%     P _t\cA_t^{\rT}
%    +
%    \cB_t\cB_t^{\rT},\\
%\label{dotR1}
%    \dot{ R }_t
%    & = &
%    \cA_t R _t
%    +
%     R _t\cA_t^{\rT}
%    +
%     P _t\cC_t^{\rT}\cC_t P _t,
%\end{eqnarray}
%with initial conditions  $ P _0 = P$ and $ R _0 = R$. The
%variational problem (\ref{JT}) is an optimal control problem for the
%deterministic dynamical system governed by
%(\ref{dotP1})--(\ref{dotR1}) in the state space $(\mS_{2n}^+)^2$.
%This is  a time-varying covariance control problem.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section{Hamilton-Jacobi-Bellman equation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%Assuming that the optimal cost function (\ref{JT}) is continuously
%differentiable with respect to $(T,P,R) \in \mR_+ \x \mS_{2n}^+ \x
%\mS_{2n}^+$ in the sense of Frechet, it satisfies the
%Hamilton-Jacobi-Bellman (HJB) equation
%\begin{eqnarray}
%\nonumber
%    \d_TJ_T
%    & = &
%    \min_{a,b,c}
%    (
%        \Tr
%        (
%            \cC
%            (P+\theta R)
%            \cC^{\rT}
%        )
%        \\
%\nonumber
%        & & \, \quad+ \,
%        \Tr
%        (
%            \Phi
%            (
%                \cA P+ P \cA^{\rT}+ \cB\cB^{\rT}
%            )
%        )\\
%\label{HJB}
%        & & \, \quad+\,
%        \Tr
%        (
%            \Psi
%            (
%                \cA R+ R\cA^{\rT}+ P\cC^{\rT}\cC P
%            )
%        )
%    ),
%\end{eqnarray}
%where
%\begin{equation}
%\label{PhiPsi}
%    \Phi
%    :=
%    \d_P J_T,
%    \qquad
%    \Psi:= \d_R J_T
%\end{equation}
%are the Frechet derivatives of the function $J_T$. The function
%under minimization on the right-hand side of (\ref{HJB}) is affine
%in $a$. More precisely, its dependence on $a$ is described by
%$2\Tr((\Phi P + \Psi R)_{11}a^{\rT})$. Therefore, for the minimum
%with respect to $a$ to be finite, it is necessary that
%\begin{equation}
%\label{PR11}
%    \mho_{11} = 0,
%\end{equation}
%where
%\begin{equation}
%\label{mho}
%    \mho
%    :=
%    \Phi P + \Psi R.
%\end{equation}
% Indeed, if (\ref{PR11}) is not satisfied, then the
%minimization with respect to  $a$ in (\ref{HJB}) yields  $-\infty$.
%Moreover, under the smoothness of $J_T$ in $P$ and $R$, the relation
%(\ref{PR11}) also follows from the invariance
%\begin{equation}
%\label{invariance}
%    J_T(D_S(P), D_S(R)) = J_T(P,R),
%\end{equation}
%which holds for any nonsingular matrix $S$ of order $n$, describing
%the controller state transformation (\ref{abc_transform}). Here,
%$D_S$ is the associated transformation of the controllability
%Gramian of the closed-loop system:
%\begin{equation}
%\label{DS}
%    D_S(P)
%    =
%    \left[
%        \begin{array}{cl}
%            S & 0\\
%            0 & I_n
%        \end{array}
%    \right]
%    P
%    \left[
%        \begin{array}{ll}
%            S^{\rT} & 0\\
%            0 & I_n
%        \end{array}
%    \right].
%\end{equation}
%Indeed, by differentiating the left-hand side of (\ref{invariance})
%with respect to $S$ at $S=I_n$ for fixed but otherwise arbitrary $P$
%and $R$ partitioned according to (\ref{blocks}), it follows that
%\begin{eqnarray*}
%    0
%    & = &
%    \left.\d_S J_T(D_S(P),D_S(R))\right|_{S=I_n}\\
%    & = &
%    2(
%    \Phi_{1\bullet}P_{\bullet 1}
%    +
%    \Psi_{1\bullet}R_{\bullet 1})
%     =
%    2\mho_{11}.
%\end{eqnarray*}
%
%
%
%
%%==============================================================================
%\begin{lemma}
%\label{lem:bc} Suppose the matrix $D_{21}$ is of full row rank,
%and $D_{12}$ is of full column rank. Also, let the state-space
%realization $(a,b,c)$ of a stabilizing controller,
%satisfying (\ref{full:cA})--(\ref{full:dLdc}), be minimal. Then
%\begin{eqnarray}
%\label{b}
%    b
%    & = &
%    -\Phi_{11}^{-1}
%    (
%        \mho_{12}C_2^{\rT}
%        +
%        \Phi_{12} B_1 D_{21}^{\rT}
%    )
%    (
%        D_{21}
%        D_{21}^{\rT}
%    )^{-1},\\
%\nonumber
%    c
%    & = &
%    -
%    (
%        D_{12}^{\rT}
%        D_{12}
%    )^{-1}
%    (
%        B_2^{\rT}
%        \mho_{21}
%        +
%        D_{12}^{\rT}C_1(P+\theta R + P\Psi P)_{21}
%    )\\
%\label{c}
%    & \x &
%    (P+\theta R + P\Psi P)_{11}^{-1}.
%\end{eqnarray}
%\end{lemma}
%%==============================================================================
%
%\textit{Proof of Lemma~\ref{lem:bc}:} Substitution of
%the matrices $\cB$ and $\cC$ from (\ref{full:cB}) and
%(\ref{full:cC}) into (\ref{full:dLdb}) and (\ref{full:dLdc}) brings
%these equations to the form
%\begin{eqnarray}
%%\label{full:dLdb1}
%    \mho_{12}C_2^{\rT}
%    +
%    (
%        \Phi_{11}bD_{21}
%        +
%        \Phi_{12}B_1
%    )
%    D_{21}^{\rT}
%    & = &
%    0,\\
%%\label{full:dLdc1}
%    B_2^{\rT}\mho_{21}
%    +
%    D_{12}^{\rT}
%    (
%        D_{12}c (P+\theta R +P\Psi P)_{11}
%        +
%        C_1(P+\theta R +P\Psi P)_{21}
%    )
%    & = &
%    0.\qquad
%\end{eqnarray}
%By Lemma~\ref{lem:pos_def_11}, the matrices $\Phi_{11}$ and
%$\Psi_{11}$ are nonsingular. Therefore, left multiplication of both
%sides of (\ref{full:dLdb1}) by $\Phi_{11}^{-1}$ and right
%multiplication by $(D_{21} D_{21}^{\rT})^{-1}$ yields (\ref{b}).
%Similarly, right multiplication of both sides of (\ref{full:dLdc1})
%by $\Psi_{11}^{-1}$ and left multiplication by $(D_{12}^{\rT}
%D_{12})^{-1}$ yields (\ref{c}). \endproof
%%==============================================================================
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section{Solution via Pontryagin's minimum principle}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%The Hamiltonian function associated with the covariance control
%problem is
%\begin{eqnarray}
%\nonumber
%    H(\underbrace{P,R}_{\rm state},\underbrace{a,b,c}_{\rm control},\underbrace{\Phi,\Psi}_{\rm costate})
%    & = &
%        \Tr
%        (
%            \cC
%            (P+\theta R)
%            \cC^{\rT}
%        )
%        \\
%\nonumber
%        & + &
%        \Tr
%        (
%            \Phi
%            (
%                \cA P+ P \cA^{\rT}+ \cB\cB^{\rT}
%            )
%        )\\
%%\label{H}
%        & + &
%        \Tr
%        (
%            \Psi
%            (
%                \cA R+ R\cA^{\rT}+ P\cC^{\rT}\cC P
%            )
%        ).
%\end{eqnarray}
%Here, the dependence on the controller state-space realization
%matrices $a$, $b$ and $c$ is accomplished through the matrices
%$\cA$, $\cB$ and $\cC$. Application of Pontryagin's minimum
%principle \cite{?} yields the time-reversed equations for the
%costate $(\Phi, \Psi)$ of the covariance system
%\begin{eqnarray}
%\label{dotPhi}
%    -\dot{\Phi}
%    & = &
%    \d_P H
%    =
%    \cA^{\rT}\Phi + \Phi \cA +
%    \cC^{\rT}\cC
%    +\Psi P \cC^{\rT}\cC + \cC^{\rT}\cC P\Psi,\\
%\label{dotPsi}
%    -\dot{\Psi}
%    & = &
%    \d_R H
%    =
%    \cA^{\rT}\Psi + \Psi \cA + \theta \cC^{\rT}\cC.
%\end{eqnarray}
%
%%$$
%%    \min \bQ_{\theta,T}
%%    =
%%    \Tr(\Phi_0 P_0)
%%    +
%%    HT
%%    -
%%    \int_{0}^{T}
%%        \Tr(P\d_P H + R\d_R H)
%%    \rd t
%%$$
%
%
%
%
%\newpage
%
%With the matrix
%$\cA$ assumed Hurwitz, let $ P $ and $ Q $ denote the
%controllability and observability Gramians of the closed-loop system
%(\ref{F}) satisfying the algebraic Lyapunov equations
%\begin{equation}
%\label{cPQ}
%    \cA P
%    +
%     P \cA^{\rT}
%    +
%    \cB\cB^{\rT}
%    = 0,
%\quad
%    \cA^{\rT} Q
%    +
%     Q  \cA
%    +
%    \cC^{\rT}\cC
%    = 0.
%\end{equation}
%As before, we will utilize the matrix $H$ associated with $P$ and $Q$ by (\ref{H}).
%In the subsequent optimization, $ P $ and $ Q $ are employed as
%additional independent variables with values in the space of real
%symmetric matrices of order $2n$ which are subject to the linear
%constraints (\ref{cPQ}). By the last representation of the
%quadro-quartic functional (\ref{QQstate}), the Lagrange function,
%associated with the LQQG problem (\ref{LQQG}), is
%\begin{eqnarray}
%\nonumber
%    \Lambda(a,b,c, P , Q ,\Phi,\Psi)
%    :=
%    \Tr
%    (
%        \Phi
%        (
%            \cA P  +  P \cA^{\rT} + \cB\cB^{\rT}
%        )
%    )\\
%\nonumber
%    +
%    \Tr
%    (
%        \Psi
%        (
%            \cA^{\rT} Q  +  Q  \cA  + \cC^{\rT}\cC
%        )
%    )\\
%\label{Lagrange}
%    -
%    2\Tr
%    (
%        \cA^{\rT}
%         H
%        (I+\theta  H)
%    ).
%\end{eqnarray}
%Here, the state-space realization matrices $\cA$, $\cB$, $\cC$ of
%the closed-loop system are affine functions of the controller
%matrices $a$, $b$, $c$ described by (\ref{ABC}). Also, $\Phi$ and
%$\Psi$ are real symmetric matrices of order $2n$ which play the role
%of Lagrange multipliers associated with the linear constraints
%(\ref{cPQ}).
%
%To obtain equations for a critical point of the
%Lagrange function (\ref{Lagrange}), we compute its Frechet
%derivatives with respect to $a$, $b$, $c$, $ P $, $ Q $, $\Phi$,
%$\Psi$ and equate them to zero. Equating the Frechet derivatives of
%$\Lambda$ with respect to the Lagrange multipliers $\Phi$ and $\Psi$
%in (\ref{Lagrange}) replicates the Lyapunov equations (\ref{cPQ}):
%\begin{eqnarray}
%\label{dLdPhi}
%    \d_{\Phi} \Lambda
%    & = &
%        \cA  P
%        +
%         P  \cA^{\rT}
%    +
%    \cB\cB^{\rT} = 0,\\
%\label{dLdPsi}
%    \d_{\Psi} \Lambda
%     & = &
%    \cA^{\rT} Q
%    +
%     Q  \cA
%    +
%    \cC^{\rT}\cC = 0.
%\end{eqnarray}
%
%
%%\begin{eqnarray}
%%\label{dLda}
%%    \d_a \Lambda
%%    & = &
%%    2
%%    (
%%        \d_{\cA}\Lambda
%%    )_{11} = 0,\\
%%\label{dLdb}
%%    \d_b \Lambda
%%    & = &
%%    2
%%    \left(
%%        \d_{\cA} \Lambda
%%    \right)_{12}C_2^{\rT}
%%    +
%%    2
%%    \Phi_{1\bullet}\cB D_{21}^{\rT} = 0,\\
%%\label{dLdc}
%%    \d_c \Lambda
%%    & = &
%%    2B_2^{\rT}
%%    (
%%        \d_{\cA} \Lambda
%%    )_{21}
%%    +
%%    2D_{12}^{\rT}\cC
%%    \Psi_{\bullet 1} = 0,\\
%%%\nonumber\\
%%%\nonumber
%%%    \d_{ P } \Lambda
%%%    & = &
%%%    \cA^{\rT} (\Phi- Q -\theta Q  P  Q )
%%%    +
%%%    (\Phi- Q -\theta Q  P  Q ) \cA\\
%%%\label{dLdP}
%%%    & - &
%%%    \theta
%%%         Q (\cA  P  +  P \cA^{\rT}) Q
%%%    = 0,\\
%%%%\nonumber\\
%%%\nonumber
%%%    \d_{ Q } \Lambda
%%%    & = &
%%%    \cA(\Psi- P -\theta  P  Q  P )  +
%%%    (\Psi- P -\theta  P  Q  P )\cA^{\rT}\\
%%%\label{dLdQ}
%%%     & - &
%%%    \theta
%%%         P (\cA^{\rT}  Q  +  Q \cA) P
%%%     = 0,
%%%     \\
%%%\nonumber\\
%%%\label{dLdPhi}
%%%    \d_{\Phi} \Lambda
%%%    & = &
%%%        \cA  P
%%%        +
%%%         P  \cA^{\rT}
%%%    +
%%%    \cB\cB^{\rT} = 0,\\
%%%\label{dLdPsi}
%%%    \d_{\Psi} \Lambda
%%%    & = &
%%%    \cA^{\rT} Q
%%%    +
%%%     Q  \cA
%%%    +
%%%    \cC^{\rT}\cC = 0.
%%\end{eqnarray}
%%====================================================================
%
%%====================================================================
%
%%====================================================================
%\begin{lemma}
%\label{lem:PhiPsi} At a critical point of the Lagrange function
%$\Lambda$ from (\ref{Lagrange}), with $\cA$ Hurwitz, the Lagrange
%multipliers $\Phi$ and $\Psi$ are expressed in terms of the Gramians
%$ P $ and $ Q $ from (\ref{cPQ}) and the matrix $H$ from (\ref{H}) as
%\begin{eqnarray}
%    \Phi
%     & = &
%     Q  + \theta ( H  Q  + \wt{\Phi}),\\
%    \Psi
%     & = &
%     P  + \theta ( P  H  + \wt{\Psi}).
%\end{eqnarray}
%\end{lemma}
%%====================================================================
%
%\begin{proof}
%The Frechet derivatives of the Lagrange function $\Lambda$ from
%(\ref{Lagrange}) with respect to the Gramians $ P $ and $ Q $ are
%computed as
%\begin{eqnarray}
%\nonumber
%    \d_{ P } \Lambda
%    & = &
%    \cA^{\rT} (\Phi- Q -\theta H  Q )
%    +
%    (\Phi- Q -\theta H  Q ) \cA\\
%\label{dLdP}
%    & - &
%    \theta
%         Q (\cA  P  +  P \cA^{\rT}) Q ,\\
%%\nonumber\\
%\nonumber
%    \d_{ Q } \Lambda
%    & = &
%    \cA(\Psi- P -\theta  P  H )  +
%    (\Psi- P -\theta  P  H )\cA^{\rT}\\
%\label{dLdQ}
%     & - &
%    \theta
%         P (\cA^{\rT}  Q  +  Q \cA) P .
%\end{eqnarray}
%The assertion of the lemma is now obtained by equating these
%derivatives to zero and substituting $\cA  P  +  P \cA^{\rT} =
%-\cB\cB^{\rT}$ from (\ref{dLdPhi}) into (\ref{dLdP}) and $\cA^{\rT}
% Q  +  Q \cA = -\cC^{\rT}\cC$ from (\ref{dLdPsi}) into (\ref{dLdQ}).
%\end{proof}
%%====================================================================


\end{document}
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



%==============================================================================
\begin{lemma}
\label{lem:bilinear} Let $\alpha$ and $\beta$ be $\mR^s$-valued
random vectors,  and let $\xi$ and $\eta$ be $\mR^t$-valued random
vectors, with the four vectors being jointly Gaussian distributed
with zero mean. Then
\begin{equation}
\label{bilinear}
    \cov(
        \alpha^{\rT} \beta,
        \xi^{\rT} \eta
    )
    =
    \Tr(\cov(\alpha,\xi)\cov(\eta,\beta) + \cov(\alpha,\eta)\cov(\xi,\beta)).
\end{equation}
\end{lemma}
%==============================================================================

\textit{Proof of Lemma~\ref{lem:bilinear}:}
Applying the representation \cite{Isserlis_1918} for the mixed
 moment of Gaussian random
variables in terms of their covariances to the entries of the random
vectors $\alpha$, $\beta$, $\xi$ and $\eta$,
\begin{eqnarray*}
    \bE(
        \alpha_i \beta_i \xi_j \eta_j
    )
    & = &
    \bE(\alpha_i\beta_i)
    \bE(\xi_j\eta_j)\\
    & + &
    \bE(\alpha_i\xi_j)
    \bE(\eta_j \beta_i )\\
    & + &
    \bE(\alpha_i\eta_j)
    \bE(\xi_j\beta_i ).
\end{eqnarray*}
Therefore,
\begin{eqnarray*}
    \bE(
        \alpha^{\rT} \beta \xi^{\rT} \eta
    )
    & = &
    \sum_{1\< i\<s,\ 1\<j\< t}
    \bE(\alpha_i \beta_i \xi_j \eta_j)
     =
    \bE(\alpha^{\rT} \beta)
    \bE(\xi^{\rT} \eta)\\
    & + &
    \Tr(\cov(\alpha,\xi)\cov(\eta,\beta) +
    \cov(\alpha,\eta)\cov(\xi,\beta)),
\end{eqnarray*}
whence (\ref{bilinear}) follows.
\endproof
%==============================================================================



%==============================================================================
\begin{lemma}
\label{lem:QR} Let $\xi$ and $\eta$ be jointly Gaussian random
vectors with zero mean and values in $\mR^s$ and $\mR^t$,
respectively, and let $Q =Q^{\rT}\in \mR^{s\x s}$ and $R =
R^{\rT}\in \mR^{t\x t}$. Then
\begin{equation}
\label{QR}
    \cov(
        \xi^{\rT} Q\xi,
         \eta^{\rT}R\eta)
    )
    =
    2
    \Tr
    (Q \cov(\xi,\eta)R\cov(\eta,\xi)).
\end{equation}
\end{lemma}
%==============================================================================

\textit{Proof:} Applying the representation \cite{Isserlis_1918} for the mixed
 moment of Gaussian random
variables in terms of their covariances to the entries of the random
vectors $\xi$ and $\eta$,
\begin{eqnarray*}
    \bE(\xi_i\xi_j\eta_k\eta_{\ell})
    & = &
    \bE(\xi_i\xi_j)\bE(\eta_k\eta_{\ell})\\
    & + &
    \bE(\xi_i\eta_k)\bE(\eta_{\ell} \xi_j )\\
    & + &
    \bE(\xi_i\eta_{\ell})\bE(\eta_k \xi_j).
\end{eqnarray*}
Therefore,
\begin{eqnarray*}
    \bE(\xi^{\rT} Q \xi \eta^{\rT} R\eta)
    & = &
    \sum_{1\< i,j\<s,\ 1\<k,\ell\< t}
    q_{ij}r_{kl}
    \bE(\xi_i\xi_j\eta_k\eta_l)\\
    & = &
    \bE(\xi^{\rT} Q \xi)
    \bE(\eta^{\rT} R\eta)
    +
    2
    \Tr
    (Q \cov(\xi,\eta)R\cov(\eta,\xi)),
\end{eqnarray*}
whence (\ref{QR}) follows.
\endproof
%==============================================================================



For example, if both $Q$ and $R$ are the identity
matrices, then Lemma~\ref{lem:QR} yields
\begin{equation}
\label{cov_energy}
    \cov(|\xi|^2, |\eta|^2)
    =
    2\|\cov(\xi,\eta)\|^2.
\end{equation}
In particular,
\begin{equation}
\label{var_energy}
    \var(|\xi|^2)
    =
    2\|\cov(\xi)\|^2.
\end{equation}
Furthermore, stacking $\xi$ and $\eta$ one underneath the other in
the direct sum $\zeta:= (\xi^{\rT}, \eta^{\rT})^{\rT}$, so that
$|\zeta|^2 = |\xi|^2 + |\eta|^2$, and using (\ref{cov_energy}),
\begin{equation}
\label{var_sum_energy}
    \var(|\zeta|^2)
    =
    \var(|\xi|^2)+\var(|\eta|^2)
    +
    4\|\cov(\xi,\eta)\|^2.
\end{equation}
Thus, the variance of the squared Euclidean norm $|\cdot|^2$ of a
subvector of a zero mean Gaussian random vector is {\it
superadditive} with respect to the direct sum of vectors.

%Using higher order mixed moments, similarly to the proof of
%Lemma~\ref{lem:QR},
%\begin{equation}
%\label{E3}
%    \bE(|\xi|^6)
%    =
%    8\Tr(C^3) + 6\Tr C\, \Tr(C^2) + (\Tr C)^3,
%    \qquad
%    C:= \cov(\xi).
%\end{equation}
%Thus, writing
%$$
%    E_k
%    :=
%    \bE(|\xi|^{2k})
%$$
%for the $k$th moment of $|\xi|^2$ associated with the Gaussian
%random vector $\xi \sim \cN(0,C)$, the relations (\ref{var_energy})
%and (\ref{E3}) give
%\begin{eqnarray*}
%E_1
%    & = &
%    \Tr C,\\
%E_2
%    & = &
%    2\Tr (C^2) + E_1^2,\\
%E_3
%    &= &
%     8\Tr(C^3)  - 2E_1^3 + 3E_1 E_2.
%\end{eqnarray*}
%This chain of equalities can be infinitely continued, and an
%additional analysis shows that
%$$
%    E_k = (2k-2)!! \Tr(C^k) + P_k(E_1, \ldots, E_{k-1}).
%$$
%Here, $(2r)!! = r!2^r$, and $P_k$ is a polynomial of the preceding
%moments $E_1, \ldots, E_{k-1}$ which is   homogeneous in the sense
%that it is a linear combination of products
% $\prod_{j=1}^{k-1} E_j^{a_j}$, where $a_1, \ldots,
%a_{k-1}$ are nonnegative integers satisfying $\sum_{j=1}^{k-1} ja_j
%= k$.
%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section*{Appendix~A. Verification of the $\cH_4$-norm property}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\renewcommand{\theequation}{A.\arabic{equation}}
%\setcounter{equation}{0}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%Let $F \in \mR^{p\x m}$ be an arbitrary matrix. In addition to its
%Frobenius norm $\|F\|$, define
%\begin{equation}
%\label{F4}
%    \|F\|_4
%    :=
%    \sqrt[4]{\Tr ((FF^{\rT})^2)}.
%\end{equation}
%We will show that $\|\cdot\|_4$ is a norm on $\mR^{p\x m}$. Apart
%from the nonegativeness and positive homogeneity which are obvious
%for $\|\cdot\|_4$ as well as the identification property , the only
%property whose verification is nontrivial is convexity. For this
%purpose, we will calculate the second order Frechet derivative of
%$\|F\|_4$ at a fixed but otherwise arbitrary $F\ne 0$ and verify its
%positive semi-definiteness. The first derivative is
%$$
%    \|F\|_4'
%    =
%    \|F\|_4^{-3}
%    FF^{\rT} F.
%$$
%Differentiating again, the second order Frechet derivative of
%$\|F\|_4$ is a self-adjoint linear operator on $\mR^{p \x m}$ given
%by
%\begin{eqnarray*}
%    \|F\|_4''(G)
%    & = &
%    \|F\|_4^{-3}
%        (
%            GF^{\rT} F +
%            FG^{\rT} F +
%            FF^{\rT} G
%        )\\
%    & - &
%    3
%    \|F\|_4^{-7}
%    \Bra
%        FF^{\rT} F,
%        G
%    \Ket
%    FF^{\rT} F,
%\end{eqnarray*}
%where $\bra G, H\ket := \Tr(G^{\rT} H)$ is the Frobenius inner
%product of real matrices. Therefore,
%\begin{eqnarray*}
%    \Bra
%        G,
%        \|F\|_4''(G)
%    \Ket
%    & = &
%    \|F\|_4^{-7}\\
%    & \x &
%    \left(
%          \|F\|_4^4
%          \Tr(
%              G^{\rT}
%              (
%                  GF^{\rT} F +
%                  FG^{\rT} F +
%                  FF^{\rT} G
%              )
%          ) -
%          3
%          \Bra
%              FF^{\rT} F,
%              G
%          \Ket^2
%    \right)
%\end{eqnarray*}
%Now, assuming,  without loss of generality, that $p\< m$, let
%$$
%    F = USV^{\rT}
%$$
%be a singular value decomposition of $F$, where $U$ and $V$ are
%orthogonal matrices of order $p$ and $m$, respectively, and
%$$
%    S
%    =
%    \left[
%        \begin{array}{c|c}
%            \Sigma & 0_{p\x (m-p)}
%        \end{array}
%    \right],
%    \qquad
%    \Sigma
%    :=
%    \diag_{1\< k \< p} (\sigma_k),
%$$
%where $\sigma_1, \ldots, \sigma_p$ are the square roots of the
%eigenvalues of $FF^{\rT}=U\Sigma^2 U^{\rT}$, so that
%$$
%    \|F\|_4 = \sqrt[4]{\sum_{k=1}^{p}\sigma_k^4}
%$$


The state-space realization of $F$ is therefore
\begin{equation}
\label{realF}
    F
     :=
    \left[
        \begin{array}{cc}
            F_{11} & F_{12} \\
            F_{21} & F_{22}
        \end{array}
    \right]
    =
    \left[
        \begin{array}{l|ll}
            A & B_1 & B_2 \\
            \hline
            C_1 & \ 0 & D_{12} \\
            C_2 & D_{21} & \ 0
        \end{array}
    \right],
\end{equation}
where the blocks
\begin{equation}
\label{subsystem}
    F_{ij}
    :=
    \left[
        \begin{array}{l|l}
            A & B_j \\
            \hline
            C_i & D_{ij}
        \end{array}
    \right],
    \qquad
    1 \< i,j \< 2,
\end{equation}
are subsystems of $F$.   If the control signal $U:= KY$ is produced
at the output of a controller $K$, which is a LCTI causal system
with the observation signal $Y$ at the input, then
\begin{equation}
\label{UKY}
    U = K Y,
\end{equation}
and the resulting closed-loop system with the exterior  input $W$
and output $Z$ is obtained by a linear-fractional transformation of
the pair $(F,K)$,
\begin{equation}
\label{LFT}
    \fL(F,K)
    :=
    F_{11} + F_{12} K
    (
        I_{p_2}
        - F_{22} K
    )^{-1}
    F_{21};
\end{equation}
see Fig.~\ref{fig:blocks}.
%==============================================================================
\begin{figure}[htb]
\begin{center}
%\includegraphics[width=5cm]{Images/blocks.ps}
\end{center}
\caption{
    \label{fig:blocks}
    The closed-loop system (\ref{LFT}) with external input $W$ and
    output $Z$.
}
\end{figure}
%==============================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{State estimator}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


$$
    \wh{x}_k
    :=
    \bE(x_k \mid \cF_{k-1}^Y)
$$
for the one-step predictor of the state signal $X$. Then, by the
$\cF_k^Y$-measurability of the control $u_k$,
$$
    \wh{x}_{k+1}
    =
    \left[
        \begin{array}{c|c}
            A & B_1
        \end{array}
    \right]
    \bE
    \left(
        \left[
            \begin{array}{c}
                x_k \\
                w_k
            \end{array}
        \right]
        \mid
        \cF_k^Y
    \right)
    +
    B_2u_k
$$

%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%We now develop necessary conditions for optimality of a
%stabilizing controller $K$ for the LQQG problem (\ref{LQQG}) in the
%class (\ref{K}) with a fixed controller state dimension $\nu$. With
%the matrix $\cA$ assumed Hurwitz, we denote by $ P $ and $ Q $ the
%controllability and observability Gramians of the closed-loop system
%(\ref{F}) satisfying the algebraic Lyapunov equations
%\begin{equation}
%\label{cPQ}
%    \cA P  +  P \cA^{\rT} + \cB\cB^{\rT} = 0,
%    \qquad
%    \cA^{\rT} Q  +  Q  \cA  + \cC^{\rT}\cC = 0.
%\end{equation}
%In the subsequent optimization, however, $ P $ and $ Q $ are
%employed as additional independent variables with values in the
%space of real symmetric matrices of order $n+\nu$ which are subject
%to the linear constraints (\ref{cPQ}). Dropping the factor $1/2$ in
%the symmetrized representation of the quadro-quartic functional
%(\ref{symQQstate}), the Lagrange function is
%\begin{eqnarray}
%\nonumber
%    \Lambda(a,b,c, P , Q ,\Phi,\Psi)
%    & := &
%    \Tr
%    (
%        \cC^{\rT}\cC P
%        +
%        \cB\cB^{\rT} Q
%    )\\
%\nonumber
%    & + &
%    \!\!\!\!\theta\,
%    \Tr
%    (
%        \cC^{\rT}\cC P  Q  P
%        +
%        \cB\cB^{\rT} Q  P  Q
%    )\\
%\nonumber
%    & + &
%    \Tr
%    (
%        \Phi
%        (
%            \cA P  +  P \cA^{\rT} + \cB\cB^{\rT}
%        )
%    )\\
%    & + &
%\label{Lagrange}
%    \Tr
%    (
%        \Psi
%        (
%            \cA^{\rT} Q  +  Q  \cA  + \cC^{\rT}\cC
%        )
%    ).
%\end{eqnarray}
%Here, the state-space realization matrices $\cA$, $\cB$ and $\cC$ of
%the closed-loop system
% depend affinely on the controller matrices
%$a$, $b$ and $c$ as described by (\ref{ABC}). Also, $\Phi$ and
%$\Psi$ are real symmetric matrices of order $n+\nu$ which play the
%role of Lagrange multipliers associated with the linear constraints
%(\ref{cPQ}). To obtain equations for a critical point of the
%Lagrange function (\ref{Lagrange}), we first compute its Frechet
%derivatives:
%\begin{eqnarray}
%\label{dLda}
%    \d \Lambda/\d a
%    & = &
%    2(\Phi P  +  Q  \Psi)_{22},\\
%\label{dLdb}
%    \d \Lambda/\d b
%    & = &
%    2(\Phi P  +  Q  \Psi)_{21}C_2^{\rT}
%    +
%    2(\Phi +  Q +\theta Q  P  Q )_{2\bullet}\cB D_{21}^{\rT},\\
%\label{dLdc}
%    \d \Lambda/\d c
%    & = &
%    2B_2^{\rT}(\Phi P  +  Q  \Psi)_{12}
%    +
%    2D_{12}^{\rT}\cC (\Psi+ P +\theta P  Q  P )_{\bullet 2},\\
%\nonumber\\
%\label{dLdP}
%    \d \Lambda/\d  P
%    & = &
%    \cA^{\rT} \Phi
%    +
%    \Phi \cA +
%    \cC^{\rT}\cC
%    \  +
%    \theta
%    (
%         Q \cB\cB^{\rT} Q
%        +
%        \cC^{\rT}\cC P  Q
%        +
%         Q  P \cC^{\rT}\cC
%    ),\\
%\label{dLdQ}
%    \d \Lambda/\d  Q
%    & = &
%    \cA\Psi  +
%    \Psi\cA^{\rT} +
%    \cB\cB^{\rT}
%    +
%    \theta
%    (
%         P \cC^{\rT}\cC P
%        +
%        \cB\cB^{\rT} Q  P
%        +
%         P  Q \cB\cB^{\rT}
%    ),\\
%\nonumber\\
%\label{dLdPhi}
%    \d \Lambda/\d \Phi
%    & = &
%        \cA  P
%        +
%         P  \cA^{\rT}
%    +
%    \cB\cB^{\rT},\\
%\label{dLdPsi}
%    \d \Lambda/\d \Psi
%    & = &
%    \cA^{\rT} Q
%    +
%     Q  \cA
%    +
%    \cC^{\rT}\cC.
%\end{eqnarray}
%Here, we have used the Frechet derivatives of the Lagrange function
%$\Lambda$ with respect to the matrices $\cA$, $\cB$ and $\cC$:
%\begin{eqnarray}
%\label{dLdcA}
%    \d \Lambda/\d \cA
%    & = &
%    2(\Phi P  +  Q  \Psi),\\
%\label{dLdcB}
%    \d \Lambda/\d \cB
%    & = &
%    2(\Phi+  Q +\theta Q  P  Q )\cB,\\
%\label{dLdcC}
%    \d \Lambda/\d \cC
%    & = &
%    2\cC (\Psi+ P +\theta P  Q  P ).
%\end{eqnarray}
%Also, the Gramians $ P $ and $ Q $, and related matrices $\Omega$
%(including $\Phi$ and $\Psi$), inherit the four block structure of
%the matrix $\cA$ in (\ref{ABC}). The sizes and the numbering of the
%blocks are as follows
%$$
%    \Omega
%    :=
%    \begin{array}{cc}
%    \!{}_{\leftarrow\, n\, \rightarrow}\  {}_{\leftarrow\, \nu\,\rightarrow} &\\
%        \left[
%            \begin{array}{ll}
%                \Omega_{11} & \Omega_{12}\\
%                \Omega_{21} & \Omega_{22}
%            \end{array}
%        \right]
%    &\!\!\!\!\!\!
%        \begin{array}{c}
%            \updownarrow\!{}^n\\
%            \updownarrow\!{}_{\nu}
%        \end{array}\\
%        {}
%    \end{array}
%    =
%    \begin{array}{cc}
%    \!{}_{\leftarrow\, n\, \rightarrow}\  {}_{\leftarrow\, \nu\,\rightarrow} &\\
%        \left[
%            \begin{array}{ll}
%                \Omega_{\bullet 1} & \Omega_{\bullet 2}
%            \end{array}
%        \right]
%    &\!\!\!\!\!\!
%        \begin{array}{c}
%            \updownarrow\!{}^{n+\nu}
%        \end{array}\\
%        {}
%    \end{array}
%    =
%    \begin{array}{cc}
%    \!{}_{\leftarrow n+\nu \rightarrow}\\
%        \left[
%            \begin{array}{ll}
%                \Omega_{1\bullet}\\
%                \Omega_{2\bullet}
%            \end{array}
%        \right]
%    &\!\!\!\!\!\!
%        \begin{array}{c}
%            \updownarrow\!{}^n\\
%            \updownarrow\!{}_{\nu}
%        \end{array}\\
%        {}
%    \end{array}.
%$$
%The derivatives of $\Lambda$ with respect to $\Phi$ and $\Psi$
%replicate the left-hand sides of the corresponding Lyapunov
%equations (\ref{cPQ}).
% These derivatives all vanish at a critical point of
%the Lagrange function. By equating the Frechet derivatives
%(\ref{dLdP}) and (\ref{dLdQ}) to zero, it follows that
%\begin{eqnarray}
%    \Phi
%    & = &
%     Q  + \theta  Q  P  Q  + \wt{\Phi},\\
%    \Psi
%    & = &
%     P  + \theta  P  Q  P  + \wt{\Psi},
%\end{eqnarray}
%where $\wt{\Phi}$ and $\wt{\Psi}$ are positive-semidefinite real
%symmetric  matrices satisfying the Lyapunov equations
%\begin{eqnarray}
%    \cA^{\rT}\wt{\Phi} + \wt{\Phi}\cA + 2\theta  Q \cB\cB^{\rT}  Q
%    & = & 0,\\
%    \cA\wt{\Psi} + \wt{\Psi}\cA^{\rT} + 2\theta  P \cC^{\rT}\cC  P
%    & = &
%    0.
%\end{eqnarray}
%These representations follow from the lemma below which further
%develops the Gramian anticommutation idea of Remark~2.
%
%%====================================================================
%\begin{lemma}
%\label{lem:QPQ_PQP} Let $P$ and $Q$ are the controllability and
%observability Gramians of an asymptotically stable system $(A,B,C)$.
%Then the matrices $QPQ$ and $PQP$ satisfy the Lyapunov equations
%\begin{eqnarray}
%\label{QPQ}
%    A^{\rT} QPQ + QP QA
%    & = &
%    QBB^{\rT}Q
%    -C^{\rT}C
%    PQ
%    -
%    QPC^{\rT}C,\\
%\label{PQP}
%    APQP  + PQP A^{\rT}
%    & = &
%    PC^{\rT}CP
%    -BB^{\rT}QP
%    -
%    PQBB^{\rT}.
%\end{eqnarray}
%
%\end{lemma}
%%====================================================================
%\paragraph*{Proof of Lemma~\ref{lem:QPQ_PQP}.} By applying the
%anticommutation form (\ref{comm2}) of the second Lyapunov equation
%for the observability Gramian $Q$ from (\ref{PQ}), it follows that
%\begin{eqnarray*}
%    \underbrace{A^{\rT} Q}PQ + QP\underbrace{QA}
%    & = &
%    -(QA + C^{\rT}C)
%    PQ
%    -
%    QP
%    (A^{\rT}Q + C^{\rT}C)\\
%    & = &
%    -
%    Q\underbrace{(AP+PA^{\rT})}Q
%    -C^{\rT}C
%    PQ
%    -
%    QPC^{\rT}C\\
%    & = &
%    QBB^{\rT}Q
%    -C^{\rT}C
%    PQ
%    -
%    QPC^{\rT}C,
%\end{eqnarray*}
%which proves (\ref{QPQ}).  Here, we have also used the first
%Lyapunov equation from (\ref{PQ}) for the controllability Gramian
%$P$. The other identity (\ref{PQP}) follows from (\ref{QPQ}) by the
%duality argument.
%\endproof
%%====================================================================
%$$
%    \cA
%    =
%    \left[
%        \begin{array}{cc}
%            A & 0\\
%            0 & 0
%        \end{array}
%    \right]
%    +
%    \left[
%        \begin{array}{c}
%            0\\
%            I
%        \end{array}
%    \right]
%    a
%    \left[
%        \begin{array}{cc}
%            0 & I
%        \end{array}
%    \right]
%    +
%    \left[
%        \begin{array}{c}
%            0\\
%            I
%        \end{array}
%    \right]
%    b
%    \left[
%        \begin{array}{cc}
%            C_2 & 0
%        \end{array}
%    \right]
%    +
%    \left[
%        \begin{array}{c}
%            B_2\\
%            0
%        \end{array}
%    \right]
%    c
%    \left[
%        \begin{array}{cc}
%            0 & I
%        \end{array}
%    \right]
%$$
\begin{equation}
\label{ABC}
    \cA
    :=
          \left[
              \begin{array}{cc}
                  a &  bC_2\\
                  B_2 c & A
              \end{array}
          \right],
    \qquad
    \cB
    :=
          \left[
              \begin{array}{cc}
                  bD_{21}\\
                  B_1
              \end{array}
          \right],
    \qquad
    \cC
    :=
          \left[
              \begin{array}{cc}
                  D_{12} c & C_1
              \end{array}
          \right].
\end{equation}


$$
    \cA
    =
    \left[
        \begin{array}{cc}
            A & 0\\
            0 & 0
        \end{array}
    \right]
    +
    \left[
        \begin{array}{c}
            0\\
            I
        \end{array}
    \right]
    a
    \left[
        \begin{array}{cc}
            0 & I
        \end{array}
    \right]
    +
    \left[
        \begin{array}{c}
            0\\
            I
        \end{array}
    \right]
    b
    \left[
        \begin{array}{cc}
            C_2 & 0
        \end{array}
    \right]
    +
    \left[
        \begin{array}{c}
            B_2\\
            0
        \end{array}
    \right]
    c
    \left[
        \begin{array}{cc}
            0 & I
        \end{array}
    \right]
$$

\begin{eqnarray}
\frac{\d \cA}{\d a}(\alpha)
    & = &
    \left[
        \begin{array}{c}
            0\\
            I
        \end{array}
    \right]
    \alpha
    \left[
        \begin{array}{cc}
            0 & I
        \end{array}
    \right],\\
\frac{\d \cA}{\d b}(\beta)
    & = &
    \left[
        \begin{array}{c}
            0\\
            I
        \end{array}
    \right]
    \beta
    \left[
        \begin{array}{cc}
            C_2 & 0
        \end{array}
    \right],\\
\frac{\d \cA}{\d c}(\gamma)
    & = &
    \left[
        \begin{array}{c}
            B_2\\
            0
        \end{array}
    \right]
    \gamma
    \left[
        \begin{array}{cc}
            0 & I
        \end{array}
    \right]
\end{eqnarray}


\begin{eqnarray}
\label{dLdcA}
    \d \Lambda/\d \cA
    & = &
    2(\Phi P  +  Q  \Psi -  Q  P  (I+\theta  Q  P ) ),\\
\label{dLdcB}
    \d \Lambda/\d \cB
    & = &
    2\Phi \cB,\\
\label{dLdcC}
    \d \Lambda/\d \cC
    & = &
    2\cC \Psi.
\end{eqnarray}
%More precisely, $G_{t|s}$ is a
%$\mR^{2n\x 2n}$-valued function satisfying the PDE
%$$
%    \d_tG_{t|s}=\cA_tG_{t|s},
%    \qquad
%    s\< t,
%$$
%with the boundary condition $G_{s|s}=I_{2n}$.


$$
    \d_t\var(\cE_t)
    =
    4
    \int_{0}^{t}
    \|\cov(z_s,z_t)\|^2
    \rd s
$$
For any intermediate time $0<t<T$, let
\begin{equation}
\label{EE}
    \cE_t^-
    :=
    \int_{0}^{t}
    |z_r|^2
    \rd r,
    \qquad
    \cE_t^+
    :=
    \int_{t}^{T}
    |z_s|^2
    \rd s,
\end{equation}
denote the past and future output energies of the system.

\begin{eqnarray}
\nonumber
             \left[
                 \begin{array}{cc}
                     \rho  &  \sigma \\
                     \tau & 0
                 \end{array}
             \right]
             & = &
             \left[
                 \begin{array}{c}
                     I_r\\
                     0_{p\x r}
                 \end{array}
             \right]
             \rho
             \left[
                 \begin{array}{cc}
                     I_r & 0_{r\x m}
                 \end{array}
             \right]
             +
             \left[
                 \begin{array}{c}
                     I_r\\
                     0_{p\x r}
                 \end{array}
             \right]
             \sigma
             \left[
                 \begin{array}{cc}
                     0_{m\x r} & I_m
                 \end{array}
             \right]\\
 \label{natbij}
             & + &
             \left[
                 \begin{array}{c}
                     0_{r\x p}\\
                     I_p
                 \end{array}
             \right]
             \tau
             \left[
                 \begin{array}{cc}
                     I_p & 0_{p\x r}
                 \end{array}
             \right].
\end{eqnarray}

%The natural bijection between the spaces $\bGamma_{r,m,p}$ and
%$\mR^{r\x r}\x \mR^{r\x m}\x \mR^{p\x r}$ is described by the linear
%map
%$$
%                 \begin{bmatrix}
%                     \rho  &  \sigma \\
%                     \tau & 0
%                 \end{bmatrix}
%             =
%                 \begin{bmatrix}
%                     I\\
%                     0
%                 \end{bmatrix}
%             \rho
%                 \begin{bmatrix}
%                     I & 0
%                 \end{bmatrix}
%             +
%                 \begin{bmatrix}
%                     I\\
%                     0
%                 \end{bmatrix}
%             \sigma
%                 \begin{bmatrix}
%                     0 & I
%                 \end{bmatrix}
%             +
%                 \begin{bmatrix}
%                     0\\
%                     I
%                 \end{bmatrix}
%             \tau
%                 \begin{bmatrix}
%                     I & 0
%                 \end{bmatrix}.
%$$


\begin{equation}
\label{blocks}
    \cA
    :=
    \begin{array}{cc}
    \!{}_{\leftarrow\, n\, \rightarrow}\  {}_{\leftarrow\, n\,\rightarrow} &\\
        \left[
            \begin{array}{ll}
                \cA_{11} & \cA_{12}\\
                \cA_{21} & \cA_{22}
            \end{array}
        \right]
    &\!\!\!\!\!\!
        \begin{array}{c}
            \updownarrow\!{}^n\\
            \updownarrow\!{}_n
        \end{array}\\
        {}
    \end{array}
    =
    \begin{array}{cc}
    \!{}_{\leftarrow\, n\, \rightarrow}\  {}_{\leftarrow\, n\,\rightarrow} &\\
        \left[
            \begin{array}{ll}
                \cA_{\bullet 1} & \cA_{\bullet 2}
            \end{array}
        \right]
    &\!\!\!\!\!\!
        \begin{array}{c}
            \updownarrow\!{}^{2n}
        \end{array}\\
        {}
    \end{array}
    =
    \begin{array}{cc}
    \!{}_{\leftarrow 2n \rightarrow}\\
        \left[
            \begin{array}{ll}
                \cA_{1\bullet}\\
                \cA_{2\bullet}
            \end{array}
        \right]
    &\!\!\!\!\!\!
        \begin{array}{c}
            \updownarrow\!{}^n\\
            \updownarrow\!{}_n
        \end{array}\\
        {}
    \end{array}.
\end{equation}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\subsection{Auxiliary noise shaping filter}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\renewcommand{\theequation}{C.\arabic{equation}}
%\setcounter{equation}{0}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%Let $\chi$ be a function of a complex variable, analytic in a
%neighbourhood of the origin where it vanishes: $\chi(0)=0$. Also, suppose $\chi$ is real-valued and strictly increasing on the
%interval $[0,\lambda]$ where
%$$
%    \lambda
%    :=
%    \sup
%    \{
%        \mu>0:\
%        \chi\
%        {\rm is\ analytic\ in\ a\ neighbourhood\ of}\
%        [0,\mu]
%    \}.
%$$
%We will refer to $\chi$ as the {\it gauge function} with analyticity
%radius $\lambda$. The associated functional
%\begin{equation}
%\label{Jchi}
%    \bJ_{\chi}(F)
%    :=
%    \frac{1}{2\pi}
%    \int_{-\infty}^{+\infty}
%    \Tr
%    \chi
%    (
%        \wh{F}^*
%        \wh{F}
%    )
%    \rd \omega
%\end{equation}
%(with $2\pi$ introduced for convenience) can be interpreted as a
%performance index of an asymptotically stable system $F$ whose
%$\cH_{\infty}$-norm does not exceed the analyticity radius of
%$\chi$:
%\begin{equation}
%\label{threshold}
%    \|F\|_{\infty}
%    <
%    \lambda.
%\end{equation}
%For example, the performance criterion, which is used in the
%risk-sensitive control \cite{?} and minimum entropy
%$\cH_{\infty}$-control \cite{MG_1991} , is a particular case of
%(\ref{Jchi}) with the gauge function
%\begin{equation}
%\label{RSchi}
%    \chi(v)
%    =
%    -\ln(1-\theta v),
%\end{equation}
%where $\theta>0$ is a parameter. Similarly, the quadro-quartic
%functional (\ref{QQ}) is also a special case of (\ref{Jchi}) with
%\begin{equation}
%\label{QQchi}
%    \chi(v)= v + \theta v^2/2.
%\end{equation}
%Now, let $\mF$ be a finite-dimensional manifold of asymptotically
%stable systems $F$. In particular, such is the set of closed-loop
%systems $F$ obtained for a given finite-dimensional plant $P$ with
%stabilizing controllers $K$ of a fixed state dimension.
%
%%==============================================================================
%\begin{lemma}
%\label{lem:filter} Let a  system $F\in \mF$ satisfy
%(\ref{threshold}). Suppose $G$ is an asymptotically stable system
%satisfying the factorization
%\begin{equation}
%\label{GGchi}
%    \wh{G}\wh{G}^* = \chi\,'(\wh{F}^*\wh{F}).
%\end{equation}
%Then $F$ is a critical point of the functional $\bJ_{\chi}(F)$ from
%(\ref{Jchi}) if and only if it is a critical point of the weighted
%$\cH_2$-norm $
%    \|FG\|_2
%$, where $G$ is fixed and plays the role of a noise shaping filter.
%\end{lemma}
%%==============================================================================
%
%\textit{Proof of Lemma~\ref{lem:filter}:} Let
%$
%    \chi(v)=
%    \sum_{k=1}^{+\infty}c_k v^k
%$ be the Taylor series expansion of the gauge function which shares
%the analyticity radius $\lambda$ with its derivative $
%    \chi\,'(v)
%    =
%    \sum_{k=1}^{+\infty} c_kkv^{k-1}
%$. By Lemma~\ref{lem:H2kvariation}, the first variation of the
%functional $\bJ_{\chi}(F)$ is computed as
%\begin{eqnarray}
%\nonumber
%    \delta \bJ_{\chi}(F)
%    & = &
%%    \sum_{k=1}^{+\infty}
%%    c_k
%%    \delta
%%        \|F\|_{2k}^{2k}
%%    =
%    \frac{1}{\pi}
%    \int_{-\infty}^{\infty}
%    \Tr
%    (
%        \wh{F}
%        \sum_{k=1}^{+\infty}
%        c_k
%        k
%        \Phi^{k-1}
%        \wh{\delta F}^*
%    )
%    \rd \omega\\
%\nonumber
%    & = &
%    \frac{1}{\pi}
%    \int_{-\infty}^{\infty}
%    \Tr
%    (
%        \wh{F}
%        \chi\,'(\Phi)
%        \wh{\delta F}^*
%    )
%    \rd \omega\\
%\label{deltaJ}
%    & = &
%    \frac{1}{\pi}
%    \int_{-\infty}^{\infty}
%    \Tr
%    (
%        \wh{F}
%        \wh{G}\wh{G}^*
%        \wh{\delta F}^*
%    )
%    \rd \omega
%    =
%    \delta (\|FG\|_2^2).
%\end{eqnarray}
%Here, use has been made of the factorization (\ref{GGchi}). The
%variation in (\ref{deltaJ}) is considered for the fixed noise
%shaping filter $G$, though the latter is specified by $F$. Since
%$\delta\bJ_{\chi}(F)$ and $\delta(\|FG\|_2^2)$ vanish
%simultaneously, hence, the assertion of Lemma~\ref{lem:filter}
%follows.
%\endproof



%%====================================================================
%\begin{lemma}
%\label{lem:QPQ_PQP} Let $P$ and $Q$ be the controllability and
%observability Gramians of an asymptotically stable system $(A,B,C)$.
%Then the matrices $QPQ$ and $PQP$ satisfy the Lyapunov equations
%\begin{eqnarray}
%\label{QPQ}
%    A^{\rT} QPQ + QP QA
%    & = &
%    QBB^{\rT}Q
%    -C^{\rT}C
%    PQ
%    -
%    QPC^{\rT}C,\\
%\label{PQP}
%    APQP  + PQP A^{\rT}
%    & = &
%    PC^{\rT}CP
%    -BB^{\rT}QP
%    -
%    PQBB^{\rT}.
%\end{eqnarray}
%
%\end{lemma}
%%====================================================================
%
%\textit{Proof of Lemma~\ref{lem:QPQ_PQP}:} By applying the
%anticommutation form (\ref{comm2}) of the second Lyapunov equation
%for the observability Gramian $Q$ from (\ref{PQ}), it follows that
%\begin{eqnarray*}
%    \underbrace{A^{\rT} Q}PQ + QP\underbrace{QA}
%    & = &
%    -(QA + C^{\rT}C)
%    PQ
%    -
%    QP
%    (A^{\rT}Q + C^{\rT}C)\\
%    & = &
%    -
%    Q\underbrace{(AP+PA^{\rT})}Q
%    -C^{\rT}C
%    PQ
%    -
%    QPC^{\rT}C\\
%    & = &
%    QBB^{\rT}Q
%    -C^{\rT}C
%    PQ
%    -
%    QPC^{\rT}C,
%\end{eqnarray*}
%which proves (\ref{QPQ}).  Here, we have also used the first
%Lyapunov equation from (\ref{PQ}) for the controllability Gramian
%$P$. The other identity (\ref{PQP}) follows from (\ref{QPQ}) by the
%duality argument.
%\endproof
%%====================================================================
%==============================================================================
%%==============================================================================
%\begin{figure}[htb]
%\begin{center}
%\includegraphics[width=3cm]{Images/zfw.ps}
%\end{center}
%\vskip-3mm
%\caption{
%    \label{fig:zfw}
%    A LTI system $F$ with input $W$ and
%    output $Z$.
%}
%\end{figure}
% 
1208.3815.extracted_bib
\begin{thebibliography}{99}
%--------------------------------------------------------------------
\bibitem{AM_1990}
B.D.O.Anderson, and J.B.Moore, {\it Optimal Control: Linear
Quadratic Methods}, Prentice Hall, Englewood Cliffs, 1990.
%%--------------------------------------------------------------------
%\bibitem{Avram_1988}
%F.Avram, On bilinear forms in Gaussian random variables and
%Toeplitz matrices, {\it Probab. Theory Relat. Fields}, vol. 79,
%1988, pp. 37--45.
%%--------------------------------------------------------------------
%\bibitem{CR_2007}
%C.D.Charalambous and  F.Rezaei, ``Stochastic uncertain systems
%subject to relative entropy constraints: induced norms and
%monotonicity properties of minimax games", {\it IEEE Trans. Autom.
%Control}, vol. 52, no. 4, pp. 647--663, Apr. 2007.
%--------------------------------------------------------------------
%\bibitem{DE_1997}
%P.Dupuis, and R.S.Ellis, {\it A Weak Convergence Approach to the
%Theory of Large Deviations}, Wiley, 1997.
%--------------------------------------------------------------------
\bibitem{DJP_2000}
P.Dupuis, M.R.James, and I.R.Petersen, Robust properties of
risk-sensitive control, {\it Math. Control Signals Systems}, vol.
13, 2000, pp. 318--332.
%--------------------------------------------------------------------
\bibitem{Frisch_1995}
U.Frisch, {\it Turbulence: The Legacy of A.N.Kolmogorov}, Cambridge
University Press, New York, 1995.
%--------------------------------------------------------------------
\bibitem{Ginovian_1994}
M.S.Ginovian, On Toeplitz type quadratic functionals of stationary
Gaussian processes, {\it Probab. Theory Relat. Fields}, vol. 100, 1994,
pp. 395--406.
%====================================================================
%\bibitem{Girsanov_1960}
%I.V.Girsanov, ``On transforming a certain class of stochastic
%processes by absolutely continuous substitution of measures'', {\it
%Theor. Probability Appl.}, vol. 5, no. 3, pp. 285--301, 1960.
%--------------------------------------------------------------------
%\bibitem{Gray_2008}
%R.M.Gray, {\it Entropy and Information Theory}, Springer-Verlag, New
%York,  2008.
%--------------------------------------------------------------------
\bibitem{GS_1958}
U.Grenander, and G.Szeg\H{o}, {\it Toeplitz Forms and Their
Applications}, University of California Press, 1958.
%--------------------------------------------------------------------
\bibitem{HJ_2007}
 R.A.Horn, and C.R.Johnson, {\it Matrix Analysis}, Cambridge
University Press, New York, 2007.
%--------------------------------------------------------------------
\bibitem{Isserlis_1918}
L.Isserlis, On a formula for the product-moment coefficient of any
order of a normal frequency distribution in any number of
variables, {\it Biometrika}, vol. 12, no. 1/2, 1918, pp. 134--139.
%--------------------------------------------------------------------
%\bibitem{Kailath_1969}
%T.Kailath, �A General Likelihood-Ratio Formula for Random Signals in
%Gaussian Noise�, {\it IEEE Trans. Information Theory}, vol. 15, no.
%3, pp. 350--361, 1969.
%--------------------------------------------------------------------
\bibitem{Kimura_1997}
H.Kimura, \textit{Chain-Scattering Approach to $\cH^{\infty}$-Control},
Birkh\"{a}user, Boston, 1997.
%--------------------------------------------------------------------
%\bibitem{Kolmogorov}
%A.N.Kolmogoroff, ``${\rm \ddot{U}}$ber die analytischen Methoden in
%der Wahrscheinlichkeitsrechnung'', {\it Ann. Math.},vol. 104, pp.
%415--458, 1931.
%--------------------------------------------------------------------
\bibitem{Magnus_1978}
J.R.Magnus,  The moments of products of quadratic forms in normal
variables, {\it Statistica Neerlandica}, vol. 32, 1978, pp. 201�210.
%==============================================================================
\bibitem{M_1988}
J.R.Magnus, \textit{Linear Structures}, Oxford University Press, New York, 1988.
%--------------------------------------------------------------------
%\bibitem{MIS_1978}
%N.N.Moiseev, Yu.P.Ivanilov, and E.M.Stolyarova, {\it Optimization
%Methods}, Nauka, Moscow, 1978.
%--------------------------------------------------------------------
\bibitem{MG_1991}
D.Mustafa, and K.Glover, {\it Minimum Entropy $\cH_{\infty}$
Control}, LNCIS, vol. 146, Springer, 1991.
%--------------------------------------------------------------------
%\bibitem{UP_2001}
% V.A.Ugrinovskii and I.R.Petersen, ``Minimax LQG
%control of stochastic partially observed uncertain systems", {\it
%SIAM J. Control Optim.}, vol. 40, no. 4, pp. 1189--1226, 2001.
%--------------------------------------------------------------------
%\bibitem{SVK_1994}
%A.V.Semyonov, I.G.Vladimirov, and A.P.Kurdyukov, ``Stochastic
%approach to $H_{\infty}$-optimization'',  {\it Proc. 33rd CDC},
%Florida, USA, December 14--16, vol.~3, pp.~2249--2250, 1994.
%%--------------------------------------------------------------------
%\bibitem{Shiryaev}
%A.N.Shiryaev, {\it Probability}, Springer, New York, 1995.
%--------------------------------------------------------------------
\bibitem{Simon_2005}
B.Simon, \textit{Trace Ideals and Their Applications}, 2nd Ed.,
American Mathematical Society, Providence, RI, 2005.
%%--------------------------------------------------------------------
%\bibitem{SIG_1998}
%R.E.Skelton, T.Iwasaki, and K.M.Grigoriadis, {\it A Unified
%Algebraic Approach to Linear Control Design}, Taylor \& Francis,
%London, 1998.
%--------------------------------------------------------------------
%\bibitem{VKS_1995}
%I.G.Vladimirov, A.P.Kurdjukov, and A.V.Semyonov, ``Anisotropy of
%signals and entropy of linear time invariant systems'',  {\it
%Doklady Akad. Nauk}, vol.~342, no.~3, pp.~583--585, 1995 (in
%Russian).
%--------------------------------------------------------------------
\bibitem{VKS_1996}
I.G.Vladimirov, A.P.Kurdyukov, and A.V.Semyonov, ``On computing the
anisotropic norm of linear discrete-time-invariant systems'',  {\it
Proceedings of the 13th IFAC World Congress}, San-Francisco,
California, USA, June 30--July 5, Vol. G, 1996, pp.~179--184.
%--------------------------------------------------------------------
%\bibitem{VDK_2006}
%I.G.Vladimirov, P.Diamond, and P.Kloeden, ``Anisotropy-based robust
%performance analysis of finite horizon linear discrete time varying
%systems''. {\it Automation and Remote Control}, no.~8,
%pp.~1265--1282, 2006.
%--------------------------------------------------------------------
%\bibitem{DVKS_2001}
%P.Diamond, I.G.Vladimirov, A.P.Kurdyukov, and A.V.Semyonov,
%Ani\-so\-tro\-py-based performance analysis of linear discrete
%time invariant control systems, {\it Int. J. Contr.}, vol.~74,
%no.~1, 2001, pp.~28--42.
%--------------------------------------------------------------------
%\bibitem{Willems_1972}
%J.C.Willems, Dissipative dynamical systems. Part I: general
%theory, {\it Arch. Rational Mech. Anal.}, vol. 45, no. 22, 1972, pp.
%321--351.
%%====================================================================
%\bibitem{Arnold}
%L.Arnold, G.C.Papanicolaou, and V.Wihstutz, ``Asymptotic analysis of
%the Lyapunov exponent and rotation number of the random oscillator
%and applications'', SIAM J. Appl. Math.  {\bf 46} (1986), no.~3,
%427--450.
%%====================================================================
%\bibitem{Bird}
%R.B.Bird, W.E.Stewart, and E.N.Lightfoot, {\it Transport Phenomena},
%2nd Ed., Wiley, New York, 2002.
%%====================================================================
%\bibitem{Fedoryuk}
%M.V.Fedoryuk, {\it Asymptotic Analysis}, Springer, Berlin, 1992.
%%====================================================================
%\bibitem{Harrison}
% P.Harrison,  {\it Quantum Wells, Wires and Dots: Theoretical and Computational
% Physics}, Wiley, New  York, 2000.
%%====================================================================
%\bibitem{Jacod}
%J.Jacod and A.N.Shiryaev, {\it Limit Theorems for Stochastic
%Processes},  , 2nd Ed., Springer, Berlin, 2003.
%%====================================================================
%\bibitem{Karatzas}
%I.Karatzas and S.E.Shreve, {\it Brownian Motion and Stochastic
%Calculus}, Springer-Verlag, New York, 1988.
%%====================================================================
%\bibitem{Khasminskii}
%R.Z.Khasminskii, ``Necessary and sufficient conditions for
%asymptotic stability of linear stochastic systems'', Teor.
%Veroyatnost. i Primenen. {\bf 12} (1967), 167--172 (in Russian)
%(English translation: Theor. Probability Appl. {\bf 12} (1967),
%144--147).
%%====================================================================
%\bibitem{Khasminskii_1969}
%R.Z.Khasminskii,  {\it Stability of Systems of Differential
%Equations with Random Parametric Excitation}, Nauka, Moscow, 1969
%(in Russian).
%
%%====================================================================
%\bibitem{Kolkunov}
%V.A.Kolkunov, ``Elimination of divergences in quantum mechanics'',
%Teor. i Matem. Fiz. {\bf 3} (1970), no.~1, 72--77 (in Russian).
%%====================================================================
%\bibitem{Landau}
%L.D.Landau and E.M.Lifshits, {\it Quantum Mechanics:
%Non-relativistic Theory}, 4th Ed.,  Nauka, Moscow, 1989  (in
%Russian).
%%%====================================================================
%%\bibitem{Laparo}
%%K.A.Laparo and G.L.Blankenship, ``Almost sure instability of a class
%%of linear stochastic systems with jump process coefficients'',
%%Lecture Notes in Mathematics 1186. Lyapunov exponents (L.Arnold and
%%V.Wihstutz, eds.),  Springer, Berlin, 1986, pp.~161--190.
%%====================================================================
%\bibitem{Oseledec}
%V.I.Oseledec, ``A multiplicative ergodic theorem. Characteristic
%Lyapunov numbers for dynamical systems'', Trans. Moscow Math. Soc.
%{\bf 19} (1969), 197--232 (Translated from Russian: Trudy Moskov.
%Mat. Obsch. {\bf 19} (1968), 179--210).
%%====================================================================
%\bibitem{Pastur}
%L.Pastur and A.Figotin, {\it Spectra of Random and Almost-Periodic
%Operators},  Springer-Verlag, Berlin, 1992.
%%====================================================================
%\bibitem{Pastur_1974}
%L.Pastur, and E.P.Fel'dman, ``Wave transmittance for a thick layer
%of a randomly inhomogeneous medium'',
% Soviet Physics JETP, {\bf 40} (1974), 241--243.
%%====================================================================
%\bibitem{Pinsky}
%M.A.Pinsky, ``Instability of the harmonic oscillator with small
%noise'', SIAM J. Appl. Math. {\bf 46} (1986), no.~3, 451--463.
%%====================================================================
%\bibitem{Razavy}
%M.Razavy, {\it Quantum Theory of Tunneling}, World Scientific,
%Singapore, 2003.
%%====================================================================
%\bibitem{Rogers}
%L.C.G.Rogers and D.Williams, {\it Diffusions, Markov Processes and
%Martingales}, Cambridge University Press, Cambridge, 2000.
%%====================================================================
%\bibitem{Roy}
%D.K.Roy, {\it Quantum Mechanical Tunnelling and Its Applications},
%World Scientific, Singapore, 1986.
%%====================================================================
%\bibitem{Sakurai}
%J.J.Sakurai, {\it Modern Quantum Mechanics}, Addison-Wesley,
%Reading, Massachusetts, 1994.
%%====================================================================
%\bibitem{Shiryaev}
%A.N.Shiryaev, {\it Probability}, Springer, New York, 1996.
%%====================================================================
%\bibitem{Vladimirov}
%V.S.Vladimirov, {\it Equations of Mathematical Physics}, Mir,
%Moscow, 1984.
%%====================================================================
%\bibitem{Wihstutz}
%V.Wihstutz,  ``Parameter dependence of the Lyapunov exponent for
%linear stochastic systems. A survey'', Lecture Notes in Mathematics
%1186. Lyapunov exponents (L.Arnold and V.Wihstutz, eds.), Springer,
%Berlin, 1986, pp.~200--215.
%%====================================================================
\end{thebibliography}