# Causal, Stochastic MPC for Wave Energy Converters

Connor H. Ligeikis and Jeffrey T. Scruggs The first author was supported by an
NSF Graduate Research Fellowship. This funding is gratefully acknowledged.
Views expressed in this paper are those of the authors and do not necessarily
reflect those of the National Science Foundation.C. Ligeikis and J. Scruggs
are with the Department of Civil & Environmental Engineering, University of
Michigan, Ann Arbor, MI, 48109. Phone: 734-764-1812, email:
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We implement a causal model predictive control (MPC) strategy to maximize
power generation from a wave energy converter (WEC) system, for which the
power take-off (PTO) systems have both hard stroke (i.e., displacement) limits
and force ratings. The approach models the WEC dynamics in discrete-time, in a
manner that exactly preserves energy-flow quantities, and assumes a stationary
stochastic disturbance model for the incident wave force. The control
objective is to maximize the expected power generation in stationarity, while
accounting for parasitic losses in the power train. PTO stroke measurements
are assumed to be available for real-time feedback, as well as the free-
surface elevation of the waves at a designated location relative to the WEC,
and the open-loop dynamics of the WEC are assumed to be linear and time-
invariant. Mean-square stability of the MPC algorithm is proven. The
methodology is illustrated in a simulation example pertaining to a heaving
cylindrical buoy.

Energy Systems, Model Predictive Control, Passivity

##  1 Introduction

Wave energy converters have been researched extensively for many decades [1].
The vast array of technologies proposed to harness the resource is diverse and
defies easy generalization [2, 3, 4]. However, most WEC devices are comprised
of a mechanical system, either floating at the ocean surface or mounted to the
ocean floor, which is dynamically-excited by propagating waves [5]. This
mechanical motion is coupled to one or more PTOs (which can be either
electromechanical or hydraulic) which extract power from the structural
dynamics [6, 7, 8, 9]. The PTOs are interfaced with a power train that
transmits extracted power to a localized storage system (such as a electrical
supercapacitor, a flywheel, or a hydraulic accumulator) which serves as a
buffer, allowing for smoothed power to be transmitted to a utility grid [10,
11]. To maximize power production in irregular (i.e., stochastic) waves, there
is a tangible advantage to the control of the power extracted by the PTOs,
based on real-time feedback measurements of the WEC system’s dynamic response
[12, 13, 14, 15]. The synthesis of optimal controllers for WEC systems is a
topic that continues to be an area of active research.

Early studies on the optimal control of WECs focused primarily on linear
feedback laws. Beginning with the seminal papers by Falnes [16] and Evans
[17], it was recognized that in monochromatic waves, power is maximized by
imposing an effective impedance relationship between the PTO velocity vector,
and the opposing force or torque vector, with the optimal impedance matrix
equal to the Hermitian adjoint of the driving point impedance of the WEC
system. This result is equivalent to the impedance matching technique used in
antenna arrays and other electromagnetic technologies. Impedance matching is
of limited practicality for at least three distinct reasons. The first is
that, when extended to the case of stochastic waves, the optimal impedance-
matching feedback law is anti-causal, and therefore can only be implemented if
an accurate forecast is available for the incident wave forces on the WEC
system [18, 19]. Secondly, the technique presumes that the PTOs have no
constraints. This is especially problematic because the dynamic response of
the controlled WEC system often deviates further from equilibrium than the
uncontrolled system [20, 21, 22]. As such, it is important to account for
constraints on the PTOs, particularly their limits on displacement, and their
force ratings. Thirdly, impedance matching does not extend easily to cases in
which the WEC system exhibits nonlinear dynamics [23].

For all these reasons, the last decade has seen a rally around the use of MPC
techniques to maximize the energy generated by WEC systems [24, 21, 25, 22,
26, 27, 28, 29]. MPC is especially well-suited to address the second problem
discussed above, concerning the accommodation of PTO constraints.
Additionally, many nonlinear MPC techniques exist which can readily be applied
to WEC systems with nonlinear plant dynamics [30, 31, 32, 33]. However, most
of the MPC techniques implemented for WEC systems do not directly address the
first problem, related to causality. Indeed, many MPC techniques for WECs
assume an accurate real-time prediction exists over a long receding horizon,
for the incident wave force. Predicated on this assumption, MPC is implemented
at each discrete time step by optimizing the control force/torque trajectory
for the PTOs over this receding horizon. The first component of this optimized
trajectory is then implemented by the controller. At the next time step, the
receding horizon for the wave force forecast is advanced, the PTO input
trajectory is re-optimized, and the process is thus repeated ad infinitum.

To implement the MPC technique described above, it is necessary to somehow
forecast the incident wave force, and many techniques for this have been
proposed in the literature. Some obtain the forecast by assuming that
measurements of free-surface elevations are available at a distance up-wave
from the WEC [34]. However, most approaches presume that the incident wave
force on the WEC can be measured [35], and predict the future force
trajectories from past data using one of several algorithms, including fully-
empirical curve-fitting algorithms as well as model-based algorithms such as
the extended Kalman filter [36, 37, 38, 39, 40].

The primary purpose of this paper is to provide an alternative to the above
prediction-based techniques, and to establish that WEC MPC algorithms need not
include an explicit wave force forecast at all, to achieve near-optimal
performance. Specific contributions of the paper are:

  1. a.

We illustrate an entirely causal technique for implementing MPC in WEC
applications, which only requires that the power spectral density of the sea
state be known.

  2. b.

We illustrate that if free surface elevations are measured and available for
feedback, these measurements may be systematically incorporated into the MPC
algorithm.

  3. c.

We illustrate how parasitic losses in the PTO and power train can be
systematically incorporated into the MPC optimization objective.

  4. d.

We prove that the proposed MPC algorithm is stable in the sense of Lyapunov
(for free response), in the bounded-input bounded-state sense (for transient
disturbances), and in the mean-square sense (for stationary stochastic
disturbances).

The scope of the paper is limited to WEC systems with linear plant dynamics,
but many of the techniques we discuss may be extended to the nonlinear case.

The paper is organized as follows. Section II establishes the discrete-time
stochastic modeling framework. Section III formalizes the WEC feedback control
problem, and establishes the MPC control framework. Section IV focuses on the
trajectory optimization problem for the MPC controller, illustrates how this
algorithm can be implemented as a convex optimization, and provides the
stability results. Section V provides a numerical example of the
implementation of the MPC algorithm for a heaving, cylindrical WEC. Finally,
Section VI draws some conclusions.

###  1.1 Notation and terminology

Sets are denoted in blackboard font, e.g., $\mathbb{R}$, $\mathbb{C}$,
$\mathbb{Z}$, and so on. We notate $\mathbb{R}_{\geqslant 0}$ and
$\mathbb{Z}_{\geqslant 0}$ as the sets of all nonnegative reals and integers,
respectively. The notation $\mathbb{C}_{\leqslant 1}$ (and $\mathbb{C}_{<1}$)
denote the complex numbers with moduli less than (strictly less than) $1$.
Similar definitions hold for the opposite inequalities. For a matrix
$Q\in\mathbb{R}^{n\times m}$ or $\mathbb{C}^{n\times m}$, $Q^{T}$ and $Q^{H}$
are the transpose and Hermitian transpose, respectively. For a vector
$q\in\mathbb{R}^{n}$ or $\mathbb{C}^{n}$ and $p\in\mathbb{Z}_{\geqslant 1}$,
the norm $\|q\|_{p}=[\sum_{i=1}^{n}|q_{i}|^{p}]^{1/p}$ and the norm
$\|q\|_{\infty}=\max_{i}|q_{i}|$. For a matrix $Q\in\mathbb{R}^{n\times m}$ or
$\mathbb{C}^{n\times m}$ and $p\in\mathbb{Z}_{\geqslant 1}\cup\\{\infty\\}$,
the norm $\|Q\|_{p}$ is the induced norm, i.e.,
$\sup_{\|q\|_{p}=1}\|Qq\|_{p}$. When no subscript is specified for a norm,
i.e., $\|\cdot\|$, $p=2$ is assumed. For a matrix $Q=Q^{H}\geqslant 0$ and a
compatible vector $q$, we denote $\|q\|_{Q}=(q^{H}Qq)^{1/2}$.

For $p$ finite, the set $\mathbb{L}_{p}$ is the set of all Lebesgue-integrable
functions $f(x)$ of a real variable $x$, such that
$\|f\|_{\mathbb{L}_{p}}\triangleq(\int_{-\infty}^{\infty}\|f(x)\|^{p}dx)^{1/p}<\infty$.
The set $\mathbb{L}_{\infty}$ is the set of all Lebesgue-integrable functions
$f(x)$ such that
$\|f\|_{\mathbb{L}_{\infty}}\triangleq\sup_{x\in\mathbb{R}}\|f(x)\|<\infty$.
For a continuous-time function $f(t)$ we denote its Fourier transform as
$\hat{f}(j\omega)$. For a discrete-time function $f_{k}$, we denote its
$z$-transform as $\bar{f}(z)$. The sets $\mathbb{H}_{2}$ and
$\mathbb{H}_{\infty}$ are the Hardy spaces, comprised of complex functions
$\bar{f}(z)$ which are analytic for all $z\in\mathbb{C}_{>1}$, and which
satisfy an appropriate norm on the boundary. For $\bar{f}\in\mathbb{H}_{2}$,
we requires that
$\|\bar{f}\|_{\mathbb{H}_{2}}\triangleq(\tfrac{1}{2\pi}\int_{-\pi}^{\pi}\|\bar{f}(e^{j\Omega}\|^{2}d\Omega)^{1/2}<\infty$.
For $\bar{f}\in\mathbb{H}_{\infty}$ we require that
$\|\bar{f}\|_{\mathbb{H}_{\infty}}\triangleq\sup_{\Omega\in[-\pi,\pi]}\|\bar{f}(e^{j\Omega})\|^{2}<\infty$.

Stochastic processes and sequences, and other random variables, are denoted in
bold. For a random variable $\bm{x}$, a particular realization is denoted by
the same character in italics, i.e., $x$. For a random variable $\bm{x}$ and
information $\mathbb{I}$, we denote $\mathcal{E}\\{\bm{x}|\mathbb{I}\\}$ as
the conditional expectation of $\bm{x}$. For two random variables $\bm{x}$ and
$\bm{y}$, the expression $\mathcal{E}\\{\bm{y}|\bm{x}\\}$ is $\bm{y}$
conditioned on $\bm{x}$, and should be interpreted as a function of random
variable $\bm{x}$. With reference to a stochastic sequence
$\\{\bm{x}_{k}:k\in\mathbb{Z}\\}$ we refer to a particular realization of a
component or sub-sequence as _data_. The abbreviations $\operatorname{wp1}$
and $\operatorname{a.s.}$ stand for “with probability $1$” and “almost sure.”

##  2 System Modeling

###  2.1 Continuous-time model

For a WEC system with $n_{p}$ PTOs, let $\bm{d}(t)\in\mathbb{R}^{n_{p}}$ be
the vector of PTO displacements, and let
$\bm{v}(t)\triangleq\tfrac{d}{dt}\bm{d}(t)$ be the resultant velocity vector.
Further, let $\bm{u}(t)\in\mathbb{R}^{n_{p}}$ be the vector of colocated
forces (or torques) associated with the PTOs. We assume that $\bm{u}(t)$ can
be made to track a command with high bandwidth, and therefore can be treated
as an input that can be controlled directly. Let
$\bm{f}(t)\in\mathbb{R}^{n_{f}}$ be the vector of forces on the various
mechanical degrees of freedom of the WEC system, due to the incident waves.
Then we presume linear time-invariant (LTI) mappings $Z_{fd}$ and $Z_{ud}$,
such that

$\bm{d}(t)=\left(Z_{fd}\bm{f}\right)(t)+\left(Z_{ud}\bm{u}\right)(t).$ (1)

These mappings are uniquely characterized by their frequency response
functions $\hat{Z}_{fd}(j\omega)$ and $\hat{Z}_{ud}(j\omega)$, which are
obtained from the hydrodynamic analysis of the WEC. For simple WEC shapes,
they can be obtained via analytical series solutions to the partial
differential equations characterizing the fluid-structure interaction [41].
For more realistic WEC shapes, they can be obtained by via finite-element
techniques [42].

Let $\bm{a}(t)\in\mathbb{R}$ be the free surface elevation at some fixed
location in the ocean, not necessarily colocated with the location of the WEC
system. Then we assume $\bm{a}$ is a stationary stochastic process with known
power spectral density (PSD) $S_{a}(\omega)$ and propagatory direction, with
the convention that

$\mathcal{E}\left\\{\bm{a}(t+\tau)\bm{a}(t)\right\\}=\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{j\omega\tau}S_{a}(\omega)d\omega$
(2)

We assume linear wave theory, and consequently, that there exists a LTI
mapping $Z_{af}$ such that

$\bm{f}(t)=\left(Z_{af}\bm{a}\right)(t)$ (3)

where $Z_{af}$ is uniquely characterized by its frequency-response function
$\hat{Z}_{af}(j\omega)$, which must also be obtained through hydrodynamic
analysis. Letting $\bm{d}_{f}(t)\triangleq(Z_{fd}\bm{f})(t)$, we have that
$\bm{d}_{f}(t)$ is a stationary stochastic process as well. Defining

$\bm{y}_{f}(t)\triangleq\begin{bmatrix}\bm{d}_{f}^{T}(t)&\bm{a}(t)\end{bmatrix}^{T}$
(4)

we have that the joint spectrum of $\bm{y}_{f}$ is

$S_{y_{f}}(\omega)=\begin{bmatrix}\hat{Z}_{fd}(j\omega)\hat{Z}_{af}(j\omega)\\\
1\end{bmatrix}S_{a}(\omega)\begin{bmatrix}\hat{Z}_{fd}(j\omega)\hat{Z}_{af}(j\omega)\\\
1\end{bmatrix}^{H}.$ (5)

The stochastic system model is then, equivalently,

$\begin{bmatrix}\bm{d}(t)\\\
\bm{a}(t)\end{bmatrix}=\bm{y}_{f}(t)+\begin{bmatrix}I\\\
0\end{bmatrix}(Z_{ud}\bm{u})(t)$ (6)

Recalling that $\bm{v}(t)=\tfrac{d}{dt}\bm{d}(t)$, we have that

$\bm{v}(t)=\tfrac{d}{dt}\bm{d}_{f}(t)+\left(Z_{uv}\bm{u}\right)(t)$ (7)

where mapping $Z_{uv}$ is uniquely characterized by its frequency response
functions, as

$\displaystyle\hat{Z}_{uv}(j\omega)=$ $\displaystyle
j\omega\hat{Z}_{ud}(j\omega)$ (8)

and where we note that the spectrum for $\tfrac{d}{dt}\bm{d}_{f}$ is
$S_{v_{f}}(\omega)=\omega^{2}\hat{Z}_{fd}(j\omega)\hat{Z}_{af}(j\omega)S_{a}(\omega)\hat{Z}_{fd}^{H}(j\omega)\hat{Z}_{af}^{H}(j\omega)$.

###### Assumption 1

We assume the following:

  1. a.

$\\{\hat{Z}_{ud},\hat{Z}_{uv}\\}\subset\mathbb{L}_{\infty}$ and
$\\{S_{y_{f}},S_{v_{f}}\\}\subset\mathbb{L}_{\infty}\cap\mathbb{L}_{1}$.

  2. b.

The limit
$\hat{Z}_{uv}^{\infty}\triangleq\lim_{\omega\rightarrow\infty}\hat{Z}_{uv}(j\omega)$
exists and is real.

  3. c.

Let $\tilde{Z}_{uv}(t)$ be the impulse response of $Z_{uv}$, i.e.,

$\tilde{Z}_{uv}(t)=\hat{Z}_{uv}^{\infty}\delta(t)+\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{j\omega
t}\left(\hat{Z}_{uv}(j\omega)-\hat{Z}_{uv}^{\infty}\right)d\omega$ (9)

with a similar definition for $\tilde{Z}_{ud}(t)$. We assume
$\\{\tilde{Z}_{uv},\tilde{Z}_{ud}\\}\subset\mathbb{L}_{1}$.

  4. d.

$Z_{uv}$ (and therefore $Z_{ud}$) are causal mappings, i.e.,
$\tilde{Z}_{uv}(t)=\tilde{Z}_{ud}(t)=0$ for all $t<0$.

  5. e.

$Z_{uv}$ is a passive mapping. Equivalently, $\hat{Z}_{uv}$ is positive-real,
i.e.,

$\hat{Z}_{uv}(j\omega)+\hat{Z}_{uv}^{H}(j\omega)\geqslant 0,\ \
\forall\omega\in\mathbb{R}$ (10)

with the inequality holding strictly for all but a countably-finite subset of
$\omega\in\mathbb{R}$.

###### Remark 1

Assumptions 1a, 1b, 1c, and 1d are stated here explicitly, but are taken for
granted in most of the literature on WEC control. We note that Assumption 1c
implies that the mappings $Z_{uv}$ and $Z_{ud}$ have bounded
$\mathbb{L}_{\infty}$ gain. Because of (8), Assumption 1a necessarily implies
that $\hat{Z}_{uv}(0)=0$, and also that $\hat{Z}_{ud}$ must be strictly
proper.

###### Remark 2

If Assumption 1e is not satisfied, this implies that the WEC possesses an
internal energy source. In that case, even in the absence of a disturbance
(i.e., with $\bm{f}=0$), energy can be generated, implying the WEC is a
perpetual-motion machine. The caveat that (10) is strict for all but
countably-finite frequencies is mild. At these frequencies, it is possible to
excite the WEC with a sinusoidal force $\bm{u}(t)$ that results in zero
mechanical power injection. For almost all applications, the only (finite)
frequency where this can occur is $\omega=0$.

The power generated by the WEC at time $t$ is denoted $\bm{\varrho}(t)$, and
is assumed to be of the form

$\bm{\varrho}(t)=-\bm{u}^{T}(t)\bm{v}(t)-\bm{u}^{T}(t)R\bm{u}(t)-v_{d}^{T}\left|\bm{u}(t)\right|$
(11)

where $R\in\mathbb{R}^{n_{p}\times n_{p}}$ with $R=R^{T}>0$, and
$v_{d}\in\mathbb{R}_{\geqslant 0}^{n_{p}}$ are parameters. The first term on
the right-hand side of (11) is equal to the total power absorbed from the WEC
dynamics by the PTO at time $t$. The second and third terms are used to model
transmission losses in the power train between the PTO and the utility bus.
For an electrical power train, these two terms capture the conductive power
dissipation in the PTO and power electronics, with the second term capturing
ohmic (i.e., “$I^{2}R$”) losses, and the third term capturing diode losses.

###  2.2 Energy-preserving discretization

We assume the control input $\bm{u}(t)$ is implemented as a zero-order hold
(ZOH) signal, i.e.,

$\bm{u}(t)=\bm{u}_{k},\ \ \forall t\in[kT,(k+1)T)$ (12)

where $T$ is the sample time, and $k\in\mathbb{Z}$. We analogously refer to
continuous-time signals sampled at the transition times of the ZOH mapping via
subscripts, i.e., $\bm{d}_{k}=\bm{d}(kT)$, $\bm{a}_{k}=\bm{a}(kT)$, and so on.
Let $\bm{p}_{k}$ be the average value of $\bm{\varrho}(t)$ over the interval
$t\in[kT,(k+1)T)$. Then

$\bm{p}_{k}=-\left(\bm{u}_{k}^{T}\bm{q}_{k}+\bm{u}_{k}^{T}R\bm{u}_{k}+v_{d}^{T}\left|\bm{u}_{k}\right|\right)$
(13)

where

$\bm{q}_{k}\triangleq\tfrac{1}{T}\left(\bm{d}_{k+1}-\bm{d}_{k}\right).$ (14)

The stochastic dynamics of the resulting sampled continuous-time system can be
expressed as

$\displaystyle\bm{d}_{k}=$
$\displaystyle\left(H_{wd}\bm{w}\right)_{k}+\left(H_{ud}\bm{u}\right)_{k}$
(15) $\displaystyle\bm{a}_{k}=$ $\displaystyle\left(H_{wa}\bm{w}\right)_{k}$
(16)

where $\bm{w}_{k}\in\mathbb{R}^{n_{p}+1}$ is an independent, identically-
distributed Gaussian stochastic sequence with zero-mean and covariance
$S_{w}\geqslant 0$. Mappings $H_{ud}$, $H_{wd}$, and $H_{wa}$ are uniquely
characterized by their discrete-time frequency response functions
$\bar{H}_{ud}(e^{j\Omega})$, $\bar{H}_{wd}(e^{j\Omega})$ and
$\bar{H}_{wa}(e^{j\Omega})$. It follows that

$\displaystyle\bm{q}_{k}=$
$\displaystyle\left(H_{wq}\bm{w}\right)_{k}+\left(H_{uq}\bm{u}\right)_{k}$
(17)

where mappings $H_{wq}$ and $H_{uq}$ are uniquely characterized by their
discrete-time frequency response functions

$\displaystyle\bar{H}_{uq}(e^{j\Omega})=$
$\displaystyle\tfrac{1}{T}(e^{j\Omega}-1)\bar{H}_{ud}(e^{j\Omega})$ (18)
$\displaystyle\bar{H}_{wq}(e^{j\Omega})=$
$\displaystyle\tfrac{1}{T}(e^{j\Omega}-1)\bar{H}_{wd}(e^{j\Omega})$ (19)

In [43] it was shown that $\bar{H}_{uq}(e^{j\Omega})$, $\Omega\in[-\pi,\pi]$
is

$\bar{H}_{uq}(e^{j\Omega})=\sum\limits_{\ell=-\infty}^{\infty}\hat{Z}_{uv}(j\omega_{\ell})\operatorname{sn}^{2}(\omega_{\ell}T/2)$
(20)

where $\omega_{\ell}\triangleq\Omega/T+2\pi\ell/T$.
$\bar{H}_{ud}(e^{j\Omega})$ is then found via (18). Meanwhile, to get mappings
$\hat{H}_{wd}(e^{j\Omega})$ and $\bar{H}_{wa}(e^{j\Omega})$, we first find the
discrete-time joint PSD

$\Sigma_{y_{f}}(\Omega)=\sum\limits_{\ell=-\infty}^{\infty}\frac{1}{T}\operatorname{sn}^{2}(\omega_{\ell}T/2)S_{y_{f}}(\omega_{\ell}).$
(21)

Then, via the Spectral Factorization Theorem, we find
$\bar{H}_{wd}(e^{j\Omega})$ and $\bar{H}_{wa}(e^{j\Omega})$ as causal,
minimum-phase, discrete-time frequency-response functions satisfying

$\Sigma_{y_{f}}(\Omega)=\begin{bmatrix}\bar{H}_{wd}(e^{j\Omega})\\\
\bar{H}_{wa}(e^{j\Omega})\end{bmatrix}S_{w}\begin{bmatrix}\bar{H}_{wd}(e^{-j\Omega})\\\
\bar{H}_{wa}(e^{-j\Omega})\end{bmatrix}^{T}.$ (22)

Frequency response function $\bar{H}_{wq}(e^{j\Omega})$ is then found via
(19). We note that in (22), the factorization is not unique. Without loss of
generality, we adopt the convention of a canonical factorization [44], in
which $S_{w}$ is normalized such that

$\frac{1}{2\pi}\int_{-\pi}^{\pi}\begin{bmatrix}\bar{H}_{wd}(e^{j\Omega})\\\
\bar{H}_{wa}(e^{j\Omega})\end{bmatrix}d\Omega=I.$ (23)

Associated $z$-domain transfer functions $\bar{H}_{uq}(z)$, $\bar{H}_{wq}(z)$,
$\bar{H}_{ud}(z)$, $\bar{H}_{wd}(z)$ and $\bar{H}_{wa}(z)$ are obtained by
inverse-transforming the discrete-time frequency-response functions to get
discrete-time impulse response functions, and then $z$-transforming these. The
following theorem groups together some useful properties of these transfer
functions.

###### Theorem 1

If Assumption 1 holds then:

  1. a.

$\\{\bar{H}_{uq},\bar{H}_{ud}\\}\in\mathbb{H}_{\infty}$ and
$\\{\bar{H}_{wa},\bar{H}_{wd}\\}\subset\mathbb{H}_{2}\cap\mathbb{H}_{\infty}$.

  2. b.

$\bar{H}_{uq}$ is positive-real in discrete-time, i.e., it is analytic for all
$z\in\mathbb{C}_{>1}$ and satisfies

$\bar{H}_{uq}(z)+\bar{H}_{uq}^{H}(z)\geqslant 0,\ \ \forall
z\in\mathbb{C}_{\geqslant 1}$ (24)

with the inequality holding strictly except for $z=1$.

  3. c.

$\bar{H}_{uq}(1)=0$, and $\bar{H}_{ud}(1)+\bar{H}_{ud}^{H}(1)>0$.

###  2.3 Finite-dimensional discrete-time model

We assume that the mappings $H_{ud}$, $H_{wd}$, and $H_{wa}$ have been
approximated as finite-dimensional, LTI, discrete-time systems, using any of
several system identification techniques. Let $\bm{y}\triangleq[\bm{d}^{T}\
\bm{a}]^{T}$. Without loss of generality we assume the resulting state space
system for the mapping $\\{\bm{u},\bm{w}\\}\mapsto\bm{y}$ is in the form of an
innovations model [45], i.e.,

$H:\left\\{\begin{array}[]{rl}\bm{x}_{k+1}=&A\bm{x}_{k}+B_{u}\bm{u}_{k}+B_{w}\bm{w}_{k}\\\
\bm{y}_{k}=&C_{y}\bm{x}_{k}+\bm{w}_{k}\end{array}\right.$ (25)

and a minimal realization111Note that, as with all innovations models, the
state vector $\bm{x}_{k}$ should not be thought of as a physical state, but
rather, as the optimal state estimate for the system, conditioned on
$\\{\bm{y}_{\ell},\ell<k\\}$. Model (25) can equivalently be viewed as a
Kalman filter, with matrix $B_{w}$ as the associated Kalman gain. . For
convenience, define partitions

$C_{y}=\begin{bmatrix}C_{d}\\\ C_{a}\end{bmatrix},\quad
I=\begin{bmatrix}D_{wd}\\\ D_{wa}\end{bmatrix}$ (26)

and note that

$\displaystyle\bar{H}_{ud}(z)\approx$ $\displaystyle
C_{d}\left[zI-A\right]^{-1}B_{u}$ (27) $\displaystyle\bar{H}_{wd}(z)\approx$
$\displaystyle D_{wd}+C_{d}\left[zI-A\right]^{-1}B_{w}$ (28)
$\displaystyle\bar{H}_{wa}(z)\approx$ $\displaystyle
D_{wa}+C_{a}\left[zI-A\right]^{-1}B_{w}.$ (29)

Furthermore, $\bm{q}_{k}$ is

$\bm{q}_{k}=C_{q}\bm{x}_{k}+D_{uq}\bm{u}_{k}+D_{wq0}\bm{w}_{k}+D_{wq1}\bm{w}_{k+1}$
(30)

where

$\displaystyle C_{q}=$ $\displaystyle\tfrac{1}{T}C_{d}(A-I),$ $\displaystyle
D_{uq}=$ $\displaystyle\tfrac{1}{T}C_{d}B_{u}$ (31) $\displaystyle D_{wq0}=$
$\displaystyle\tfrac{1}{T}(C_{d}B_{w}-D_{wd}),$ $\displaystyle D_{wq1}=$
$\displaystyle\tfrac{1}{T}D_{wd}$ (32)

Consequently, we may express transfer functions $\bar{H}_{uq}$ and
$\bar{H}_{wq}$ in terms of the state space parameters as

$\displaystyle\bar{H}_{uq}(z)\approx$ $\displaystyle
D_{uq}+C_{q}[zI-A]^{-1}B_{u}$ (33) $\displaystyle\bar{H}_{wq}(z)\approx$
$\displaystyle D_{wq1}z+D_{wq0}+C_{q}[zI-A]^{-1}B_{w}.$ (34)

We presume that the finite-dimensional model has been found in a way that
preserves the properties listed in Theorem 1 for the above transfer functions.
More information on subspace-based procedures to identify the state space
parameters of $H$, while adhering to the conditions of Theorem 1, can be found
in [46] and [47]. The following theorem groups together some useful
implications of Theorem 1, which apply to the finite-dimensional model $H$ as
in (25).

###### Theorem 2

Assume the minimal finite-dimensional model (25) is such that the transfer
functions (27), (28), (29), and (33) adhere to the conditions in Theorem 1.
Then:

  1. a.

$A$ is Schur.

  2. b.

$D_{uq}+D_{uq}^{T}>0$.

  3. c.

Let $E\in\mathbb{R}^{n_{o}\times n}$ be full-row-rank with null space equal to
the unobservable subspace of $(A,C_{q})$. Then the discrete-time Riccati
equation

$W=A^{T}WA+F^{T}MF$ (35)

with

$\displaystyle M\triangleq$ $\displaystyle
R+\tfrac{1}{2}(D_{uq}+D_{uq}^{T})-B_{u}^{T}WB_{u}$ (36) $\displaystyle
F\triangleq$ $\displaystyle-M^{-1}\left[\tfrac{1}{2}C_{q}-B_{u}^{T}WA\right]$
(37)

has a solution $W=W^{T}\geqslant 0$ with $EWE^{T}>0$, $M>0$, and such that
$A+B_{u}F$ is Schur.

###  2.4 Inverse system dynamics

Consider the inverse system in which $\bm{d}$ is treated as a control input,
and $\bm{u}$ is treated as an output. Then

$\displaystyle\bm{u}_{k}=$ $\displaystyle-
D_{uq}^{-1}C_{q}\bm{x}_{k}-\tfrac{1}{T}D_{uq}^{-1}\left(C_{d}B_{w}-D_{wd}\right)\bm{w}_{k}$
$\displaystyle-\tfrac{1}{T}D_{uq}^{-1}D_{wd}\bm{w}_{k+1}+\tfrac{1}{T}D_{uq}^{-1}(\bm{d}_{k+1}-\bm{d}_{k}).$
(38)

As such,

$\bm{x}_{k+1}=\left[A-B_{u}D_{uq}^{-1}C_{q}\right]\bm{x}_{k}+\left[\tfrac{1}{T}B_{u}D_{uq}^{-1}\right]\left(\bm{d}_{k+1}-\bm{d}_{k}\right)\\\
+\tfrac{1}{T}B_{u}D_{uq}^{-1}D_{wd}(\bm{w}_{k}-\bm{w}_{k+1})\\\
+[I-\tfrac{1}{T}B_{u}D_{uq}^{-1}C_{d}]B_{w}\bm{w}_{k}.$ (39)

Changing coordinates, define $\bm{\zeta}_{k}$ such that

$\begin{bmatrix}\bm{d}_{k}\\\
\bm{\zeta}_{k}\end{bmatrix}=\begin{bmatrix}C_{d}\\\
B_{u}^{\perp}\end{bmatrix}\bm{x}_{k}+\begin{bmatrix}D_{wd}\\\
0\end{bmatrix}\bm{w}_{k}$ (40)

where $B_{u}^{\perp}$ is a full-row-rank matrix such that
$B_{u}^{\perp}B_{u}=0$ and such that $[C_{d}^{T}\ (B_{u}^{\perp})^{T}]$ is
square and invertible. (Note that this is assumed without loss of generality,
because $C_{d}B_{u}=TD_{uq}$ is nonsingular.) Then we have that

$\begin{bmatrix}C_{d}\\\
B_{u}^{\perp}\end{bmatrix}^{-1}=\begin{bmatrix}B_{u}&C_{d}^{\perp}\end{bmatrix}\begin{bmatrix}C_{d}B_{u}&0\\\
0&B_{u}^{\perp}C_{d}^{\perp}\end{bmatrix}^{-1}$ (41)

where $C_{d}^{\perp}$ is any full-column rank matrix such that
$C_{d}C_{d}^{\perp}=0$ and such that $[B_{u}\ C_{d}^{\perp}]$ is square. Using
these facts, it is straight-forward to verify that in the new coordinates,

$\displaystyle\bm{\zeta}_{k+1}=$ $\displaystyle
A_{s}\bm{\zeta}_{k}+B_{ws}\bm{w}_{k}+B_{ds}\bm{d}_{k}$ (42)
$\displaystyle\bm{x}_{k}=$ $\displaystyle
C_{s}\bm{\zeta}_{k}+D_{ws}\bm{w}_{k}+D_{ds}\bm{d}_{k}$ (43)

where

$\displaystyle A_{s}\triangleq$ $\displaystyle B_{u}^{\perp}AC_{s}$ (44)
$\displaystyle B_{ws}\triangleq$ $\displaystyle
B_{u}^{\perp}B_{w}-B_{ds}D_{wd},$ $\displaystyle B_{ds}\triangleq$
$\displaystyle B_{u}^{\perp}AD_{ds}$ (45) $\displaystyle C_{s}\triangleq$
$\displaystyle C_{d}^{\perp}\left[B_{u}^{\perp}C_{d}^{\perp}\right]^{-1}$ (46)
$\displaystyle D_{ws}\triangleq$ $\displaystyle-D_{ds}D_{wd},$ $\displaystyle
D_{ds}\triangleq$ $\displaystyle\tfrac{1}{T}B_{u}D_{uq}^{-1}$ (47)

###### Theorem 3

If the conditions of Theorems 1 and 2 hold, then $A_{s}$ is Schur.

##  3 Optimal Stochastic WEC Control

Let $\bm{Y}_{k}$ be

$\bm{Y}_{k}\triangleq\left\\{\bm{y}_{\ell}:\ell<k\right\\}$ (48)

Then we presume that at time $k$, the information available for the purposes
of control is $\bm{Y}_{k}$, i.e., we presume a control algorithm $K$ that
facilitates the mapping

$K:\bm{Y}_{k}\mapsto\bm{u}_{k}$ (49)

for each $k\in\mathbb{Z}_{\geqslant 0}$, starting from a deterministic initial
condition $\bm{x}_{0}=x_{0}$. Random variables $\bm{x}_{k}$ and $\bm{u}_{k}$
are functions of $\bm{Y}_{k}$. At time $k$ the data $\bm{Y}_{k}=Y_{k}$ is
known to the control algorithm $K$, and consequently the data
$\bm{x}_{k}=x_{k}$ and $\bm{u}_{k}=u_{k}$ are known. With each time
advancement of $k\rightarrow k+1$, innovations data $\bm{w}_{k}=w_{k}$ becomes
known, and data $\bm{x}_{k+1}=x_{k+1}$ is evaluated recursively by control
algorithm $K$, via (25).

Let $\mathbb{K}_{c}$ be the set of all strictly causal mappings $K$ as in
(49). The idealized optimal WEC control problem may be stated as the following
optimization problem:

$K=\operatorname*{sol}\left\\{\begin{array}[]{rl}\text{Given}:&A,B_{u},B_{w},C_{y},S_{w},x_{0}\\\
\text{Max}:&p_{\text{avg}}\\\ \text{Dom}:&K\in\mathbb{K}_{c}\\\
\text{Constr}:&|\bm{u}_{k}|\leqslant u_{\max}\operatorname{wp1},\ \forall
k\in\mathbb{Z}_{\geqslant 0}\\\ &|\bm{d}_{k}|\leqslant
d_{\max}\operatorname{wp1},\ \forall k\in\mathbb{Z}_{>0}\end{array}\right.$
(50)

where $u_{\max}$ and $d_{\max}$ are the vectors of force and displacement
limits for PTOs, and

$\displaystyle p_{\text{avg}}\triangleq$
$\displaystyle\lim\limits_{\tau\rightarrow\infty}\frac{1}{\tau+1}\sum\limits_{k=0}^{\tau}\mathcal{E}\\{\bm{p}_{k}\\}$
(51) $\displaystyle=$
$\displaystyle\lim\limits_{\tau\rightarrow\infty}\frac{-1}{\tau+1}\sum\limits_{k=0}^{\tau}\mathcal{E}\left\\{\bm{u}_{k}^{T}C_{q}\bm{x}_{k}+\bm{u}_{k}^{T}\tilde{R}\bm{u}_{k}+v_{d}^{T}|\bm{u}_{k}|\right\\}$
(52)

where $\tilde{R}\triangleq R+\tfrac{1}{2}(D_{uq}+D_{uq}^{T})>0$. However, we
note that optimization problem (50) is not well-posed. This is because for
each time $k$, there is a nonzero probability that there will exist no control
input that can simultaneously satisfy both the displacement and force
constraints. To rectify this, we soften the force constraint using a vector of
slack variables, denoted $\bm{b}\in\mathbb{R}^{n_{p}}_{\geqslant 0}$. We then
have the relaxed problem

$K=\operatorname*{sol}\left\\{\begin{array}[]{rl}\text{Given}:&A,B_{u},B_{w},C_{y},S_{w},x(0)\\\
\text{Max}:&p_{\text{avg}}-\lim\limits_{\tau\rightarrow\infty}\tfrac{1}{\tau+1}\sum_{k=0}^{\tau}\mathcal{E}\left\\{\mu
u_{\max}^{T}\bm{b}_{k}\right\\}\\\ \text{Dom}:&K\in\mathbb{K}_{c}\\\
\text{Constr}:&|\bm{u}_{k}|\leqslant u_{\max}+\bm{b}_{k}\operatorname{wp1},\
\forall k\in\mathbb{Z}_{\geqslant 0}\\\ &|\bm{d}_{k}|\leqslant
d_{\max}\operatorname{wp1},\ \forall k\in\mathbb{Z}_{>0}\\\
&\bm{b}_{k}\geqslant 0\operatorname{wp1},\ \forall k\in\mathbb{Z}_{\geqslant
0}\end{array}\right.$ (53)

where $\mu\in\mathbb{R}_{>0}$ is a penalty term.

The following theorem is instrumental to the formulation of this control
problem in the context of MPC, because it illustrates that the maximization of
$p_{\text{avg}}$ may be framed equivalently as the minimization of an
auxiliary positive-semidefinite objective function.

###### Theorem 4

Let $K\in\mathbb{K}_{c}$ be such that the limit

$\lim\limits_{k\rightarrow\infty}\frac{1}{k}\mathcal{E}\left\\{\|\bm{x}_{k}\|^{2}\right\\}=0$
(54)

Then

$\displaystyle p_{\text{avg}}=$
$\displaystyle\bar{p}-\lim\limits_{\tau\rightarrow\infty}\frac{1}{\tau+1}\sum\limits_{k=0}^{\tau}\mathcal{E}\left\\{\|\bm{\xi}_{k}\|_{M}^{2}+v_{d}^{T}|\bm{u}_{k}|\right\\}$
(55)

where $\bm{\xi}_{k}\triangleq\bm{u}_{k}-F\bm{x}_{k}$, matrices $M$ is as $F$
are defined in (36) and (37),
$\bar{p}\triangleq\operatorname{tr}\\{B_{w}^{T}WB_{w}S_{w}\\}$, and
$W=W^{T}\geqslant 0$ is the unique stabilizing solution to (35).

Using this theorem, we can re-express the control problem for $K$ as the
equivalent problem below:

$K=\operatorname*{sol}\left\\{\begin{array}[]{rl}\text{Given}:&A,B_{u},B_{w},C_{y},S_{w},x(0)\\\
\text{Min}:&\lim\limits_{\tau\rightarrow\infty}\frac{1}{\tau+1}\sum\limits_{k=0}^{\tau}\mathcal{E}\left\\{J(\bm{x}_{k},\bm{u}_{k},\bm{b}_{k})\right\\}\\\
\text{Dom}:&K\in\mathbb{K}_{c}\\\ \text{Constr}:&|\bm{u}_{k}|\leqslant
u_{\max}+\bm{b}_{k}\operatorname{wp1},\forall k\in\mathbb{Z}_{\geqslant 0}\\\
&|\bm{d}_{k}|\leqslant d_{\max}\operatorname{wp1},\forall
k\in\mathbb{Z}_{>0}\\\ &\bm{b}_{k}\geqslant 0\operatorname{wp1},\ \forall
k\in\mathbb{Z}_{\geqslant 0}\end{array}\right.$ (56)

where $J(x_{k},u_{k},b_{k})$ is equal to

$J(x_{k},u_{k},b_{k})\triangleq\|u_{k}-Fx_{k}\|_{M}^{2}+v_{d}^{T}|u_{k}|+\mu
u_{\max}^{T}b_{k}$ (57)

Note that because $M>0$ and $\mu>0$, $J$ is a positive-semidefinite, and is
convex.

###### Remark 3

Technically, the optimal control problem (56) is not well-posed. Satisfaction
of constraint $|\bm{d}_{k}|\leqslant d_{\max}$ must be enforced at previous
time step, $k-1$, by appropriate choice of $\bm{u}_{k-1}$. Ensuring
satisfaction of the constraint would require knowledge of innovations vectors
$\bm{w}_{k}$ and $\bm{w}_{k-1}$ and state vector $\bm{x}_{k-1}$. Due to the
fact that $K$ is constrained to be strictly causal, this knowledge is not
available to the controller until after input $\bm{u}_{k-1}$ has been applied.
As such, there is a finite (albeit very small, in practice) probability of
displacement violation, which is unavoidable. To make the optimization (56)
well-posed, the displacement constraint can be replaced with a conditional
expectation, as

$\left|\mathcal{E}\left\\{\bm{d}_{k+1}\ |\ \bm{Y}_{k}\right\\}\right|\leqslant
d_{\max}\operatorname{wp1},\forall k\in\mathbb{Z}_{\geqslant 0}$ (58)

This is what is done in the proposed MPC algorithm.

##  4 MPC Formulation

To frame the optimal stochastic WEC control problem in the context of MPC,
suppose the present time is $k$, let $h>0$ be some receding horizon length,
and presume that the data $\bm{Y}_{k}=Y_{k}$ is known. Assume that for all
time $\ell<k$, the mapping $K:\bm{Y}_{\ell}\mapsto\bm{u}_{\ell}$ is imposed,
where $K\in\mathbb{K}_{c}$. It follows that at time $k$, the data
$\bm{x}_{k}=x_{k}$ is also known.

Let the future control inputs over the receding horizon be

$\bm{U}_{k}\triangleq\\{\bm{u}_{k},...,\bm{u}_{k+h}\\}$ (59)

and let

$U_{k|k}\triangleq\\{u_{k|k},...,u_{k+h|k}\\}$ (60)

be a sequence of hypothetical, deterministic future control inputs, which are
determined from data $Y_{k}$ and therefore are probabilistically independent
of innovations $\\{\bm{w}_{k},...,\bm{w}_{k+h}\\}$. Let
$\\{x_{k|k},...,x_{k+h|k}\\}$ be the expected values of
$\\{\bm{x}_{k},...,\bm{x}_{k+h}\\}$ conditioned on data $Y_{k}$ and $U_{k|k}$,
i.e.,

$x_{m|k}\triangleq\mathcal{E}\left\\{\bm{x}_{m}\Big{|}\bm{Y}_{k}=Y_{k},\bm{U}_{k}=U_{k|k}\right\\}$
(61)

Then it follows that for $m\in\\{k,...,k+h\\}$,

$x_{m+1|k}=Ax_{m|k}+B_{u}u_{m|k}$ (62)

with initial condition $x_{k|k}=x_{k}$.

In the proposed stochastic MPC formulation of the optimal control problem, at
each time $k$ we find a deterministic input sequence $U_{k|k}$ and a
deterministic slack variable sequence

$B_{k|k}\triangleq\\{b_{k|k},...,b_{k+h,k}\\}$ (63)

which minimize an objective $\Gamma(x_{k},U_{k|k},B_{k|k})$, subject to
constraints. As such, we have that

$u_{k}=\operatorname*{sol}\limits_{u_{k|k}}\left\\{\begin{array}[]{rl}\text{Given:}&x_{k}\\\
\text{Min:}&\Gamma(x_{k},U_{k|k},B_{k|k})\\\ \text{Dom:}&U_{k|k},B_{k|k}\\\
\text{Constr:}&\Theta(x_{k},U_{k|k},B_{k|k})\leqslant 0\\\
&\Psi(x_{k},U_{k|k})=0\end{array}\right.$ (64)

for appropriately-defined constraint vector functions
$\Theta(\cdot,\cdot,\cdot)$ and $\Psi(\cdot,\cdot)$. Let the resultant
optimized trajectories be denoted by ring accents, i.e.,
$\mathring{U}_{k|k}=\\{\mathring{u}_{k|k},...,\mathring{u}_{k+h|k}\\}$ and
$\mathring{B}_{k|k}=\\{\mathring{b}_{k|k},...,\mathring{b}_{k+h|k}\\}$. Upon
solving this optimization, the MPC algorithm implements
$u_{k}=\mathring{u}_{k|k}$. Following this, data $\bm{y}_{k}=y_{k}$ is
obtained, and innovations data $\bm{w}_{k}=w_{k}$ and state data
$\bm{x}_{k+1}=x_{k+1}$ are found from (25). Time is incremented, i.e.,
$k\rightarrow k+1$, and the process is repeated ad infinitum.

###  4.1 Optimization objective

The MPC optimization objective involves the expected value of
$J(\bm{x}_{m},u_{m|k},b_{m|k})$ for $m\in\\{k,...,k+h\\}$, conditioned on data
$\bm{Y}_{k}=Y_{k}$ as well as $\bm{U}_{k}=U_{k|k}$. Because
$J(\cdot,\cdot,\cdot)$ is quadratic its first argument, the second and third
arguments are deterministic when conditioned on the data, and because
$x_{m|k}$ is an unbiased conditional estimate of $\bm{x}_{m}$, it follows that

$\mathcal{E}\left\\{\left.J(\bm{x}_{m},u_{m|k},b_{m|k})\right|\bm{Y}_{k}=Y_{k},\bm{U}_{k}=U_{k|k}\right\\}\\\
=\mathcal{E}\left\\{\|F(\bm{x}_{m}-x_{m|k})\|_{M}^{2}\right\\}+J(x_{m|k},u_{m|k},b_{m|k})$
(65)

The MPC optimization objective is the sum of this expectation over the time
horizon, plus a final-value penalty. Noting that $\bm{x}_{m}-x_{m|k}$ is zero
for $m=k$ and a linear combination of $\\{\bm{w}_{k},...,\bm{w}_{m-1}\\}$ for
$m>k$, it follows that the first term on the right-hand side of (65) is
independent of optimization variables $U_{k|k}$ and $B_{k|k}$, and can be
subtracted from the objective without affecting the solution. In so doing, the
objective to be minimized in the MPC algorithm is

$\Gamma(x_{k},U_{k|k},B_{k|k})\\\ \triangleq
J_{f}(x_{k+h+1|k})+\sum\limits_{m=k}^{k+h}J(x_{m|k},u_{m|k},b_{m|k})$ (66)

where for $m\in\\{k,...,k+h\\}$, (62) is tacitly assumed, with initial
condition $x_{k|k}=x_{k}$. Term $J_{f}(\cdot)$ is the final-value penalty,
which will be formulated in Section 4.3.

###  4.2 Constraints

Inequality constraint vector $\Theta(x_{k},U_{k|k},B_{k|k})$ is formulated to
impose the constraints on force, displacement, and the slack variables.
Specifically, we have that

$\Theta(x_{k},U_{k|k},B_{k|k})=\begin{bmatrix}|u_{k|k}|-b_{k|k}-u_{\max}\\\
\vdots\\\ |u_{k+h|k}|-b_{k+h|k}-u_{\max}\\\ |d_{k+1|k}|-d_{\max}\\\ \vdots\\\
|d_{k+h|k}|-d_{\max}\\\ -b_{k|k}\\\ \vdots\\\ -b_{k+h|k}\end{bmatrix}$ (67)

where

$\displaystyle d_{m|k}\triangleq$
$\displaystyle\mathcal{E}\left\\{\bm{d}_{m}|\bm{Y}_{k}=Y_{k},\bm{U}_{k}=U_{k|k}\right\\}=C_{d}x_{m|k}$
(68)

with (62) tacitly assumed with initial condition $x_{k|k}=x_{k}$.

Equality constraint vector $\Psi(x_{k},U_{k|k})$ plays an important role in
guaranteeing closed-loop stability of the MPC algorithm. Specifically we
impose the condition

$d_{k+h+1|k}=0.$ (69)

As such, $\Psi(x_{k},U_{k|k})=C_{d}x_{k+h+1|k}$, with (62) tacitly assumed
with initial condition $x_{k|k}=x_{k}$.

###  4.3 Final value penalty

Due to the imposition of equality constraint (69), it is known for all
$U_{k|k}$ satisfying the constraints of MPC algorithm (64), $d_{k+h+1|k}=0$.
Now, let hypothetical inputs $u_{m|k}$ beyond the time horizon, i.e., for
$m>k+h$, be those that maintain $d_{m+1|k}=0$. Using Section 2.4, this implies
post-horizon trajectories for $x_{m|k}$ and $u_{m|k}$ as

$\displaystyle x_{m|k}=$ $\displaystyle
C_{s}A_{s}^{m-(k+h+1)}B_{u}^{\perp}x_{k+h+1|k}$ (70) $\displaystyle u_{m|k}=$
$\displaystyle-D_{uq}^{-1}C_{q}x_{m|k}.$ (71)

An associated feasible post-horizon penalty term $b_{m|k}$ can then found for
all $m\geqslant k+h+1$ as the minimum value that satisfies the force
constraint, i.e.,

$b_{m|k}=\textrm{max}\\{0,|u_{m|k}|-u_{\max}\\}.$ (72)

Now, for $u_{m|k}$, and $b_{m|k}$ equal to (71) and (72) respectively,
consider the post-horizon performance measure

$\rho(x_{k+h+1|k})\triangleq\sum\limits_{m=k+h+1}^{\infty}J(x_{m|k},u_{m|k},b_{m|k})$
(73)

Ideally, we would use $\rho(\cdot)$ as the final-value penalty $J_{f}(\cdot)$
in MPC performance measure (66). However, $\rho(x_{k+h+1|k})$ is a complex
function of its argument, and in the lemma below, we show that it may be over-
bounded by a much simpler function.

###### Lemma 1

There exists a final-value penalty function $J_{f}(\cdot)$ of the form

$J_{f}(x_{k+h+1|k})=\|x_{k+h+1|k}\|_{Q_{f}}^{2}+\|C_{f}x_{k+h+1|k}\|_{1}$ (74)

where $C_{f}\in\mathbb{C}^{n_{p}\times n}$ and $Q_{f}\in\mathbb{R}^{n\times
n}$ with $Q_{f}=Q_{f}^{T}\geqslant 0$, such that the for all
$x_{k+h+1|k}\in\mathbb{R}^{n}$, the following hold:

$\displaystyle\rho(x_{k+h+1|k})\leqslant J_{f}(x_{k+h+1|k})$ (75)
$\displaystyle J_{f}(x_{k+h+2|k})+J(x_{k+h+1|k},u_{k+h+1|k},b_{k+h+1|k})$
$\displaystyle\quad\quad\quad\quad\quad\leqslant J_{f}(x_{k+h+1|k})$ (76)

where for the second inequality, $x_{m|k}$, $u_{m|k}$, and $b_{m|k}$ are
evaluated as in (70), (71), and (72) respectively.

We note that the construction of $C_{f}$ and $Q_{f}$ is given explicitly in
the proof to Lemma 1, in the appendix.

###  4.4 Stability

In this subsection we provide a stability proof for MPC algorithm (64). To be
more concise, we introduce the notation

$\displaystyle\|x\|_{1\tau}\triangleq$
$\displaystyle\sum\limits_{k=0}^{\tau}\|x_{k}\|_{1},$
$\displaystyle\|x\|_{2\tau}\triangleq$
$\displaystyle\left(\sum\limits_{k=0}^{\tau}\|x_{k}\|_{2}^{2}\right)^{1/2}$
(77)

with similar notation for norms on subsequences of $w$. Using this notation,
the theorem below distills to several common notions of stability as special
cases. For the case in which $w_{k}=0$ for all $k\geqslant 0$ and initial
condition $x_{0}\neq 0$, the theorem guarantees that $\|x\|_{2\tau}$ is
bounded from above as $\tau\rightarrow\infty$, thus guaranteeing that
$\|x_{k}\|_{2}\rightarrow 0$ as $k\rightarrow\infty$. This implies asymptotic
Lyapunov stability. For the case in which $w_{k}\neq 0$ for some $k\geqslant
0$, and $x_{0}=0$, the theorem establishes a bound on $\|x\|_{2\tau}$ which is
a continuous function of $\|w\|_{1\tau}$ and $\|w\|_{2\tau}$, which may be
interpreted as a form of bounded-input bounded-state stability. Furthermore,
if $\|\tfrac{1}{1+\tau}w\|_{1\tau}$ and $\|\tfrac{1}{1+\tau}w\|_{2\tau}$ have
finite upper bounds as $\tau\rightarrow\infty$, then the result of the theorem
guarantees a bound on
$\lim_{\tau\rightarrow\infty}\|\tfrac{1}{1+\tau}x\|_{2\tau}$. This has a
relationship to mean-square stability for the stochastic system model, as will
be shown in the corollary after the theorem.

###### Theorem 5

Let $\tau\in\mathbb{Z}_{\geqslant 0}$. Assume that for all
$k\in\\{0,...\tau\\}$, $K:Y_{k}\mapsto u_{k}$ is the MPC algorithm (64). Let
$w$ be any input sequence with $\|w_{k}\|$ bounded for all
$k\in\\{0,..,\tau\\}$ and let $x_{0}$ be any bounded initial condition. Then
there exist constants $\\{\alpha_{1}...,\alpha_{4}\\}\subset\mathbb{R}_{>0}$
such that

$\|x\|_{2\tau}^{2}\leqslant\alpha_{1}\|x_{0}\|_{1}\\\
+\alpha_{2}\|x_{0}\|_{2}^{2}+\alpha_{3}\|w\|_{1\tau}+\alpha_{4}\|w\|_{2\tau}^{2}$
(78)

For the stochastic model we consider in this paper, the above theorem may be
interpreted as follows. Consider that property (c) above holds for all
stochastic sequences $\bm{w}=w$ and all initial conditions $\bm{x}_{0}=x_{0}$
in the ensemble. Consequently, we have that

$\tfrac{1}{1+\tau}\|\bm{x}\|_{2\tau}^{2}\leqslant\tfrac{\alpha_{1}}{1+\tau}\|\bm{x}_{0}\|_{1}+\tfrac{\alpha_{2}}{1+\tau}\|\bm{x}_{0}\|_{2}^{2}+\tfrac{\alpha_{3}}{1+\tau}\|\bm{w}\|_{1\tau}+\tfrac{\alpha_{4}}{1+\tau}\|\bm{w}\|_{2\tau}^{2}$
(79)

Taking $\tau\rightarrow\infty$, we have that

$\displaystyle\lim\limits_{\tau\rightarrow\infty}\tfrac{1}{1+\tau}\|\bm{x}\|_{2\tau}^{2}\leqslant\lim\limits_{\tau\rightarrow\infty}\tfrac{1}{1+\tau}\left(\alpha_{3}\|\bm{w}\|_{1\tau}+\alpha_{4}\|\bm{w}\|_{2\tau}^{2}\right)$
(80)

But recall that each $\bm{w}_{k}$ is an independent, identically-distributed,
Gaussian random variable with zero mean and covariance $S_{w}$. Via the Strong
Law of Large Numbers, we therefore have that the following limits hold:

$\displaystyle\lim\limits_{\tau\rightarrow\infty}\tfrac{\alpha_{3}}{1+\tau}\|\bm{w}\|_{1\tau}=$
$\displaystyle\alpha_{3}\sum\limits_{i=1}^{n_{p}}\sqrt{\tfrac{2}{\pi}\hat{e}_{i}^{T}S_{w}\hat{e}_{i}},\
\ \operatorname{a.s.}$ (81)
$\displaystyle\lim\limits_{\tau\rightarrow\infty}\tfrac{\alpha_{4}}{1+\tau}\|\bm{w}\|_{2\tau}^{2}=$
$\displaystyle\alpha_{4}\operatorname{tr}\\{S_{w}\\},\ \ \operatorname{a.s.}$
(82)

We conclude that over the ensemble, there exists a bound $N$ equal to the sum
of (81) and (82), such that

$\displaystyle\lim\limits_{\tau\rightarrow\infty}\tfrac{1}{1+\tau}\|\bm{x}\|_{2\tau}^{2}\leqslant
N,\ \ \operatorname{a.s.}$ (83)

As such, we have that in a stochastic context, MPC feedback law (64) is mean-
square stable, almost-surely.

##  5 Example

To demonstrate the application of MPC algorithm (64), consider the simple WEC
system illustrated in Fig. 1. As shown, the WEC is comprised of a floating,
slack-moored, cylindrical buoy, in which a tuned vibration absorber (TVA) is
embedded. The TVA is comprised of a mass-spring-dashpot assembly, and the PTO
is situated between the mass and the buoy. The dynamics of the mass along the
buoy axis, together with the heave motion of the buoy, comprise a two degree-
of-freedom vibratory system. We presume that the mooring system restrains the
buoy motion to heave, and that the resultant dynamics can justifiably be
approximated as linear.

Figure 1: Example WEC system

The masses of the buoy and TVA were chosen such that the system is in
hydrostatic equilibrium in the configuration shown in Fig. 1. The spring of
the TVA is such that the fundamental vibratory mode of the system (including
the added mass of the displaced fluid) has a natural period of approximately
$9$s. The dashpot of the TVA is such that the fundamental mode has a fraction
of critical damping of $0.5\%$.

Regarding the PTO model, the scope of this example requires only that we
specify its parasitic loss parameters, which were chosen as
$R=1\text{kW/MN}^{2}$, and $v_{d}=1\text{kW/MN}$. (These values were deemed
reasonable if the PTO is realized as a permanent-magnet synchronous machine.)
We presume a force rating of $u_{\max}=0.5\text{MN}$ and a stroke limit of
$d_{\max}=3\text{m}$.

The free surface elevation $a(t)$ is presumed to be measured and available for
feedback. As shown in Fig. 1, it is presumed that these measurements are made
at a distance $D_{a}$ up-wave from the buoy, along the direction of wave
propagation. The specific value of $D_{a}$ will be varied in the example, with
$D_{a}=0$ implying that measurements are colocated with the buoy. As $D_{a}$
is increased from $0$, the measurements of $a(t)$ can be thought of as
providing a kind of preview of the wave loading the buoy will experience in
the near-future. However, this is only an approximate interpretation, for two
reasons. Firstly, the transfer function from the free surface elevation at the
the location of the WEC, to the corresponding wave loading force $f(t)$, is
noncausal. Consequently, even if some preview information were available for
this free-surface elevation, it would not be possible to exactly determine a
preview of $f(t)$. Secondly, ocean waves are dispersive, and consequently
waves at different temporal frequencies travel at different velocities. As
such, the time history for the free surface elevation at the buoy cannot be
viewed simply as a time-delayed version of $a(t)$ at location $D_{a}>0$.
Despite these two caveats, for the limiting case of $D_{a}\rightarrow\infty$,
knowledge of $a(\tau)$ for all $\tau<t$ implies knowledge of the future
incident wave loading $f(\tau)$ over an arbitrary but finite receding horizon
$\tau\in[t,t+t_{h}]$.

The stochastic sea state was taken to be a Pierson-Moskowitz spectrum [48],
i.e.,

$S_{a}(\omega)=c_{a}\left|\tfrac{\omega_{p}}{\omega}\right|^{5}\exp\left\\{-\tfrac{5}{4}\left(\tfrac{\omega_{p}}{\omega}\right)^{4}\right\\}$
(84)

This spectrum is traditionally parametrized by its significant wave height
$H_{s}$ and peak wave period $T_{p}$, with these parameters determining
$c_{a}$ and $\omega_{p}$ in the equation above. We have that
$T_{p}=\tfrac{2\pi}{\omega_{p}}$ and $H_{s}=4\sigma_{a}$, with
$\sigma_{a}^{2}=\tfrac{1}{2\pi}\int_{-\infty}^{\infty}S_{a}(\omega)d\omega$.
In this example, we uniformly presume that $T_{p}=9$s. We take $H_{s}=1$m
unless otherwise specified.

The sample time for the controller is uniformly taken to be $T=0.5\text{s}$.
Fig. 2 shows the discrete-time frequency response function
$\bar{H}_{uq}(e^{j\Omega})$, both for the original infinite-dimensional model
(as in Section 2.2) and the finite-dimensional approximation (as in Section
2.3). Note that the both the infinite-dimensional and finite-dimensional
models are positive-real. Analogously, Fig. 3 shows the discrete-time spectrum
$\Sigma_{y_{f}}(\Omega)$, both for the infinite-dimensional model and the
finite-dimensional approximation, showing a good match at all frequencies
where the spectrum magnitude is significant. (This spectrum corresponds to
$D_{a}=20$m.) It is straight-forward to verify from this plot that the
infinite-dimensional model adheres to Assumption 1a.

Figure 2: Frequency response function $\bar{H}_{uq}(e^{j\Omega})$. (Both
infinite-dimensional model and finite-dimensional approximation are shown, but
are indistinguishable.)  Figure 3: Spectrum $\Sigma_{y_{f}}(\Omega)$ for
infinite-dimensional model (solid) and finite-dimensional approximation
(dashed). Note that $S_{y_{f}}^{11}$ and $S_{y_{f}}^{22}$ have zero phase, and
$S_{y_{f}}^{21}$ and $S_{y_{f}}^{12}$ are complex conjugates.

MPC algorithm (64) was implemented for this model, with a time horizon of one
minute (i.e., $h=120$). The only other parameter necessary to specify is the
slack penalty factor $\mu$. Performance was found to be insensitive to this
value, so long as it is sufficiently large to prohibit violations of the force
constraint, and a value of $\mu=10$ was used for all simulations. All
simulations were performed in Matlab, and the optimizations were implemented
using CVX [49].

Fig. 4a shows the a sample path realization of the stochastic response of the
WEC for the case with $D_{a}=0$m. Note that both the force and stroke
constraints are uniformly satisfied, as desired. The generated power $p_{k}$
is also shown, and it is worth noting that the MPC algorithm requires
significant bi-directional power flow (i.e., both significantly positive and
significantly negative values). Fig. 4b shows analogous plots for the same
case but with $H_{s}=3$m. Note that although the average force magnitude is
higher than in the $H_{s}=1$m case, the force rating is only violated very
occasionally. These violations occur when it is infeasible for the MPC
algorithm to simultaneously satisfy both the stroke and force constraints.

Figure 4: Sample path realization for $a_{k}$, $d_{k}$, $u_{k}$, and $p_{k}$,
with MPC algorithm (64), for the case with $D_{a}=0$m and $H_{s}=1$m (a), and
$3$m (b)

To estimate the performance $p_{\text{avg}}$ achieved by controller, an
ensemble of $48$ sample path realizations was simulated for one hour each. The
average power generation was evaluated for each simulation, resulting in a
sample set of $48$ values. From this sample set, the mean was estimated, and
confidence intervals were assessed for this estimate. The resultant estimates
are shown in Fig. 5 for $H_{s}=1$m, and for various values of $D_{a}$. As
expected, as $D_{a}$ is increased, the generated power also increases. As
$D_{a}$ approaches $100$m, the performance asymptotically approaches a limit.

Figure 5: Estimates of $p_{\text{avg}}$ (with $90\%$ confidence intervals) vs.
$D_{a}$, with $H_{s}=1$m.

Also shown on the plot is the anticausal power generation performance. This
was obtained by optimizing $p_{\text{avg}}$ for each sample in the sample set,
assuming the entire $w$ trajectory is known _a priori_. Additionally, in this
case the optimal control trajectory was optimized all at once, rather than in
the form of an MPC algorithm. As expected, the anticausal performance exceeds
the causal performance for all values of $D_{a}$. Theoretically, the
anticausal optimal performance should be slightly higher than the asymptote of
the causal case, as $D_{a}$ is made large. This asymptotic discrepancy is due
to a small degree of sub-optimality inherent to the MPC algorithm as a
consequence of its finite horizon. However, it is worth noting that as $D_{a}$
becomes large, the anticausal and causal performance estimates are close
enough that they are within the margin of error for the simulation.

Finally, we examine the manner in which performance changes with $H_{s}$. Fig.
6 shows the estimate for $p_{\text{avg}}$ for $D_{a}=0$m, as a function of
$H_{s}$. To illustrate the features of this plot better, $p_{\text{avg}}$ is
normalized by $H_{s}^{2}$. (Note that for an unconstrained WEC with linear
dynamics and linear control, this quantity would be independent of $H_{s}$.)
Clearly, we see that, when normalized as such, the power generation
performance has a peak at approximately $H_{s}=0.5$m. Above this wave height,
the constraints on stroke and force begin to hamper the ability of the WEC to
harvest the available energy. Below this wave height, the efficiency of the
WEC power train drops significantly due to the $v_{d}^{T}|u|$ term in the
expression for the parasitic dissipation.

Figure 6: Estimates of $p_{\text{avg}}$ as a function of $H_{s}$, with
$D_{a}=0$m.

##  6 Conclusions

In this paper we have illustrated a systematic technique for synthesizing an
MPC algorithm for a stochastically-excited WEC, which adheres to force and
stroke constraints whenever possible, is provably stable, and does not require
the incident wave force to be forecast explicitly. An interesting extension to
this work would be the development of a version of the technique which can
accommodate nonlinearities in the WEC dynamics. Other items for further study
include the stability and performance robustness of the algorithm, in the
presence of uncertainty in the WEC dynamic model.

###  .1 Proof of Theorem 1

To begin, we observe that $\bar{H}_{uq}(z)$ is analytic on $\mathbb{C}_{>1}$.
Let $u^{i}_{k}=\hat{e}_{i}\delta_{k}$, where $\hat{e}_{i}$ is the Cartesian
unit vector in the $i^{th}$ direction and $\delta_{k}$ is the Kronecker delta,
and let $u^{i}(t)$ be the resultant continuous-time ZOH input, via (12). Let
$v^{i}(t)$ be the resultant continuous-time response, i.e.,
$v^{i}(t)=(Z_{uv}u^{i})(t)$. Let $q^{i}_{k}$ be the corresponding discrete-
time response. Then we note that
$\sup_{k\in\mathbb{Z}}\|q^{i}_{k}\|_{2}\leqslant\|v^{i}\|_{\mathbb{L}_{\infty}}.$
It is a standard result from input-output analysis that

$\displaystyle\|v^{i}\|_{\mathbb{L}_{\infty}}$
$\displaystyle\leqslant\|\tilde{Z}_{uv}\|_{\mathbb{L}_{1}}\|u^{i}\|_{\mathbb{L}_{\infty}}=\|\tilde{Z}_{uv}\|_{\mathbb{L}_{1}}\times
1$

which is finite due to Assumption 1c. We conclude that $\|q_{k}^{i}\|_{2}$ is
uniformly bounded for all $k\in\mathbb{Z}_{\geqslant 0}$. Now, consider that
$\bar{H}_{uq}(z)=\sum_{k=0}^{\infty}z^{-k}\begin{bmatrix}q_{k}^{1}&\cdots&q_{k}^{n_{p}}\end{bmatrix}$
which is analytic if the summation converges. Convergence is guaranteed for
$|z|>1$, due to the uniform bounds on $\|q_{k}^{i}\|_{2}$. We conclude that
$\bar{H}_{uq}(z)$ is analytic for $z\in\mathbb{C}_{>1}$. Analogous arguments
hold for $\bar{H}_{ud}(z)$. $\bar{H}_{wd}(z)$, and $\bar{H}_{wa}(z)$ are
analytic in $\mathbb{C}_{>1}$ by definition, because they are spectral
factors.

To prove claim (a), it is already known that each function is analytic for
$z\in\mathbb{C}_{>1}$, and it remains to verify that the required frequency-
domain norms hold. For $\bar{H}_{uq}$, we require that

$\|\bar{H}_{uq}\|_{\mathbb{H}_{\infty}}=\sup\limits_{\Omega\in[-\pi,\pi]}\|\bar{H}_{uq}(e^{j\Omega})\|_{2}<\infty.$

Letting $\sigma_{\ell}\triangleq\operatorname{sn}(\omega_{\ell}T/2)$, we have
that

$\displaystyle\|\bar{H}_{uq}(e^{j\Omega})\|_{2}=$
$\displaystyle\left\|\sum\limits_{\ell=-\infty}^{\infty}\hat{Z}_{uv}(j\omega_{\ell})\sigma_{\ell}^{2}\right\|_{2}\leqslant\|\hat{Z}_{uv}\|_{\mathbb{L}_{\infty}}\sum\limits_{\ell=-\infty}^{\infty}\sigma_{\ell}^{2}$

But $\sum_{\ell=-\infty}^{\infty}\sigma_{\ell}^{2}=1$ so we conclude that
$\|\bar{H}_{uq}\|_{\mathbb{H}_{\infty}}\leqslant\|\hat{Z}_{uv}\|_{\mathbb{L}_{\infty}}$.
Per Assumption 1a, we conclude that $\bar{H}_{uq}\in\mathbb{H}_{\infty}$. To
prove that $\bar{H}_{ud}\in\mathbb{H}_{\infty}$, we have that

$\displaystyle\|\bar{H}_{ud}(e^{j\Omega})\|_{2}$
$\displaystyle=\left\|\frac{1}{e^{j\Omega}-1}\sum\limits_{\ell=-\infty}^{\infty}\hat{Z}_{uv}(j\omega_{\ell})\sigma_{\ell}^{2}\right\|_{2}$
$\displaystyle=\left\|\sum\limits_{\ell=-\infty}^{\infty}\hat{Z}_{ud}(j\omega_{\ell})\left|\sigma_{\ell}\right|\right\|_{2}$
$\displaystyle\leqslant\left\|\hat{Z}_{ud}(j\omega_{0})\right\|_{2}|\sigma_{\ell}|+\sum\limits_{\begin{subarray}{c}\ell=-\infty\\\
\neq
0\end{subarray}}^{\infty}\left\|\hat{Z}_{uv}(j\omega_{\ell})\right\|_{2}\left|\frac{\sigma_{\ell}}{\omega_{\ell}}\right|$
$\displaystyle\leqslant\left\|\hat{Z}_{ud}\right\|_{\mathbb{L}_{\infty}}+\left\|\hat{Z}_{uv}\right\|_{\mathbb{L}_{\infty}}\sum\limits_{\ell=1}^{\infty}\frac{T}{(\pi\ell)^{2}}$

Due to Assumption 1a, the above bound is finite, implying that
$\|\bar{H}_{ud}\|_{\mathbb{H}_{\infty}}$ is finite.

To prove that $\bar{H}_{wa}\in\mathbb{H}_{2}$, equations (21) and (22) imply

$\displaystyle 2\pi\|\bar{H}_{wa}\|_{\mathbb{H}_{2}}^{2}$
$\displaystyle=\sum\limits_{\ell=-\infty}^{\infty}\int_{-\pi}^{\pi}\frac{1}{T}\sigma_{\ell}^{2}S_{a}(\omega_{\ell})d\Omega\leqslant\int_{-\infty}^{\infty}S_{a}(\omega)d\omega$

which is finite due to Assumption 1a. A similar argument proves the same
property for $\bar{H}_{wd}$. To prove that
$\bar{H}_{wa}\in\mathbb{H}_{\infty}$, it suffices to show that
$\|\bar{H}_{wa}(e^{j\Omega})\|_{2}$ is uniformly bounded over
$\Omega\in[-\pi,\pi]$. Again using (21) and (22),

$\displaystyle\|\bar{H}_{wa}(e^{j\Omega})\|_{2}^{2}\leqslant$
$\displaystyle\sum\limits_{\ell=-\infty}^{\infty}\frac{1}{T}\sigma_{\ell}^{2}S_{a}(\omega_{\ell})\leqslant\frac{1}{T}\|S_{a}\|_{\mathbb{L}_{\infty}}^{2}$

which is finite due to Assumption 1a. A similar argument proves the same
property for $\bar{H}_{wd}$.

For claim (b), observe that for $\Omega\in[-\pi,\pi]$,

$\bar{H}_{uq}(e^{j\Omega})+\bar{H}_{uq}^{H}(e^{j\Omega})=\sum\limits_{\ell=-\infty}^{\infty}\left[\hat{Z}_{uv}(j\omega_{\ell})+\hat{Z}_{uv}^{H}(j\omega_{\ell})\right]\sigma_{\ell}^{2}$

which is positive-semidefinite as per Assumption 1e. Because $\bar{H}_{uq}(z)$
is analytic on $\mathbb{C}_{>1}$, $\bar{H}_{uq}(z^{-1})$ is analytic on
$\mathbb{C}_{<1}$. Define

$\bar{G}_{uq}(z)\triangleq\left[I-\epsilon\bar{H}_{uq}(z^{-1})\right]\left[I+\epsilon\bar{H}_{uq}(z^{-1})\right]^{-1}$

where $0<\epsilon<\|\bar{H}_{uq}\|_{\mathbb{H}_{\infty}}^{-1}$. We then have
that (24) holds if and only if $\|\bar{G}_{uq}(z)\|\leqslant 1$ for all
$z\in\mathbb{C}_{\leqslant 1}$. But if $\bar{H}_{uq}(z^{-1})$ is analytic on
$\mathbb{C}_{<1}$ then so is $\bar{G}_{uq}(z)$, and therefore via the Maximum
Modulus Theorem, $\|\bar{G}_{uq}(z)\|$ attains is maximum on the boundary of
$\mathbb{C}_{\leqslant 1}$, i.e., for $z=e^{j\Omega}$, $\Omega\in[-\pi,\pi]$.
Because $\bar{H}_{uq}(e^{j\Omega)}+\bar{H}_{uq}(e^{-j\Omega})\geqslant 0$, it
follows that $\sup_{\Omega\in[-\pi,\pi]}\|\bar{G}_{uq}(e^{j\Omega})\|\leqslant
1$. This proves (24).

To prove the last caveat in claim (b), first consider the case where $|z|=1$,
i.e., $z=e^{j\Omega}$ for $\Omega\in[-\pi,\pi]$. Due to Assumption 1e,
evaluation of $\bar{H}_{uq}(e^{j\Omega})$ as in (20) gives that inequality
(24) must hold strictly unless $\operatorname{sn}(\omega_{\ell}T/2)=0$ for all
$\ell\in\mathbb{Z}\setminus\\{0\\}$. This only occurs if $\Omega=0$, i.e.,
$z=1$. Now, consider the case in which $z\in\mathbb{C}_{>0}$. Let
$\nu\in\mathbb{C}^{n_{p}}$ and define
$\bar{h}(z)\triangleq\nu^{H}\bar{H}_{uq}(z)\nu$. Consider the situation in
which there exists a $z_{0}\in\mathbb{C}_{>1}$ such that
$\textrm{Re}\left\\{\bar{h}(z_{0})\right\\}=0$, and define analytic function
$\bar{g}(z)\triangleq\bar{h}(z)-j\textrm{Im}\\{\bar{h}(z_{0})\\}$. Then
$\bar{g}(z_{0})=0$. Suppose it is the case that $\bar{g}(z)\neq 0$ for almost
all $z$ in an infinitesimal neighborhood of $z=z_{0}$. Let $z=z_{0}(1+\delta
e^{j\theta})$ where $\delta\in\mathbb{R}_{>0}$ is infinitesimal and
$\theta\in[-\pi,\pi]$. Because $\bar{g}(z)$ is analytic at $z=z_{0}$, and
$\delta$ is infinitesimal,

$\bar{g}(z)=\left(z-z_{0}\right)^{m}\psi(z_{0})=\left(z_{0}^{m}\delta^{m}\psi(z_{0})\right)e^{jm\theta}$

where

$\psi(z_{0})=\lim\limits_{z\rightarrow z_{0}}\frac{\bar{g}(z)}{(z-z_{0})^{m}}$

and $m\in\mathbb{Z}_{>0}$ is the lowest integer such that $\psi(z_{0})\neq 0$.
It has been proven that
$\textrm{Re}\\{\bar{g}(z)\\}=\textrm{Re}\\{\bar{h}(z)\\}\geqslant 0$ for all
$z\in\mathbb{C}_{>1}$ and therefore we require that

$z_{0}^{m}\psi(z_{0})e^{jm\theta}+z_{0}^{mH}\psi^{H}(z_{0})e^{-jm\theta}\geqslant
0,\ \ \forall\theta\in[-\pi,\pi]$ (85)

Because $z_{0}\psi(z_{0})\neq 0$, this is impossible. We therefore conclude
that if $\bar{g}(z_{0})=0$ for some $z_{0}\in\mathbb{C}_{>1}$ then it must be
the case that $\bar{g}(z)=0$ in some open neighborhood of $z=z_{0}$. But in
order for $\bar{g}(z)$ to be analytic, this implies that $\bar{g}(z)=0$ for
all $z\in\mathbb{C}_{>1}$. If this is true then $\bar{g}(z)=0$ almost
everywhere on the boundary of $\mathbb{C}_{>1}$, i.e., for $z=e^{j\Omega}$
with $\Omega\in[-\pi,\pi]$. This in turn implies that
$\textrm{Re}\\{\bar{h}(e^{j\Omega})\\}=0$ for almost all
$\Omega\in[-\pi,\pi]$, which contradicts the last caveat in assumption 1e. We
therefore conclude that $\textrm{Re}\\{\bar{h}(z)\\}>0$ everywhere in
$\mathbb{C}_{>1}$. This holds for all
$\nu\in\mathbb{C}^{n_{p}}\setminus\\{0\\}$, implying that
$\bar{H}_{uq}(z)+\bar{H}_{uq}^{H}(z)$ is nonsingular over the same domain.

The first part of claim (c) is immediate from (18), together with claim (a).
For the second part of the claim, we make a similar argument as above,
resulting in the requirement that (85) hold at $z_{0}=1$. However, since
$z_{0}=1$ is on the boundary of $\mathbb{C}_{\geqslant 1}$, rather than in the
interior, it is only necessary that condition (85) hold for
$\theta\in[-\tfrac{\pi}{2},\tfrac{\pi}{2}]$. Satisfaction of the condition
requires $m=1$ and $z_{0}\gamma(z_{0})\in\mathbb{R}_{>0}$. This proves that
$\nu^{H}\left(\bar{H}_{ud}(1)+\bar{H}_{ud}^{H}(1)\right)\nu\in\mathbb{R}_{>0}$.
Enforcing this for all $\nu\in\mathbb{C}^{n_{p}}$ implies that
$\bar{H}_{ud}(1)+\bar{H}_{ud}^{H}(1)>0$.

###  .2 Proof of Theorem 2

To prove claim (a), it is known from Theorem 1 that
$\left\\{\bar{H}_{ud},\bar{H}_{wa},\bar{H}_{wd}\right\\}\subset\mathbb{H}_{\infty}$.
Meanwhile, the mapping $\bm{u}\mapsto\bm{a}$ is zero. As such, mapping
$\\{\bm{u},\bm{w}\\}\mapsto\\{\bm{d},\bm{a}\\}$ is in $\mathbb{H}_{\infty}$.
Minimality of realization (25) then implies that $A$ is Schur.

To prove claim (b), it is a standard fact that
$D_{uq}=\tfrac{1}{2\pi}\int_{-\pi}^{\pi}\bar{H}_{uq}(e^{j\Omega})d\Omega$ and
consequently

$\frac{1}{2\pi}\int_{-\pi}^{\pi}\left[\bar{H}_{uq}(e^{j\Omega})+\bar{H}_{uq}^{H}(e^{j\Omega})\right]d\Omega=D_{uq}+D_{uq}^{T}$

The integrand is positive-definite for almost all $\Omega\in[-\pi,\pi]$ so the
integral is positive-definite.

To prove claim (c), first partition the state space to isolate the subspace
that is unobservable from $q$, i.e.,

$\displaystyle A=$ $\displaystyle\begin{bmatrix}A_{11}&0\\\
A_{21}&A_{22}\end{bmatrix},$ $\displaystyle B_{u}=$
$\displaystyle\begin{bmatrix}B_{1}\\\ B_{2}\end{bmatrix},$ $\displaystyle
C_{q}=$ $\displaystyle\begin{bmatrix}C_{1}^{T}\\\ 0\end{bmatrix}^{T}$

Presume that the theorem is known to hold for the observable subspace, i.e.,
there exists a matrix $W_{11}=W_{11}^{T}$ such that

$W_{11}=A_{11}^{T}W_{11}A_{11}+F_{1}MF_{1}^{T}$

where

$\displaystyle M=$ $\displaystyle
R+\tfrac{1}{2}\left(D_{uq}+D_{uq}^{T}\right)-B_{1}^{T}W_{11}B_{1}$
$\displaystyle F_{1}=$
$\displaystyle-M^{-1}\left[\tfrac{1}{2}C_{1}-B_{1}^{T}W_{11}A_{11}\right]$

Then the equations in the theorem are satisfied with
$W=\textrm{blockdiag}\left\\{W_{11},0\right\\}$. Henceforth we presume the
state space has been reduced by eliminating the unobservable subspace.

Now assuming the reduced system to be observable, we change coordinates to
isolate the controllable subspace, i.e.,

$\displaystyle A=$ $\displaystyle\begin{bmatrix}A_{11}&A_{12}\\\
0&A_{22}\end{bmatrix},$ $\displaystyle B_{u}=$
$\displaystyle\begin{bmatrix}B_{1}\\\ 0\end{bmatrix},$ $\displaystyle C_{q}=$
$\displaystyle\begin{bmatrix}C_{1}^{T}\\\ C_{2}^{T}\end{bmatrix}^{T}$

where we note that $A_{11}$ and $A_{22}$ are both Schur, because $A$ is.
Partitioning $W$ similarly, we have that Riccati equation (35) is equivalent
to the partitioned equations

$\displaystyle W_{11}=$ $\displaystyle A_{11}^{T}W_{11}A_{11}+F_{1}^{T}MF_{1}$
(86) $\displaystyle W_{12}=$ $\displaystyle
A_{11}^{T}W_{11}A_{12}+A_{11}^{T}W_{12}A_{22}+F_{1}^{T}MF_{2}$ (87)
$\displaystyle W_{22}=$ $\displaystyle
A_{12}^{T}W_{11}A_{12}+A_{22}^{T}W_{12}^{T}A_{12}$ (88)
$\displaystyle+A_{12}^{T}W_{12}A_{22}+A_{22}^{T}W_{22}A_{22}+F_{2}^{T}MF_{2}$
(89)

where

$\displaystyle F_{1}=$
$\displaystyle-M^{-1}\left[\tfrac{1}{2}C_{1}-B_{1}^{T}W_{11}A_{11}\right]$
(90) $\displaystyle F_{2}=$
$\displaystyle-M^{-1}\left[\tfrac{1}{2}C_{2}-B_{1}^{T}W_{11}A_{12}-B_{1}^{T}W_{12}A_{22}\right]$
(91) $\displaystyle M=$ $\displaystyle
R+\tfrac{1}{2}(D_{uq}+D_{uq}^{T})-B_{u}^{T}W_{11}B_{u}$ (92)

Noting that the system

$\left[\begin{array}[]{c|c}A_{11}&B_{1}\\\ \hline\cr
C_{1}&D_{uq}\end{array}\right]$

is a minimal realization of positive-real transfer function $\bar{H}_{uq}$,
the transfer function
$\bar{H}^{\prime}_{uq}\triangleq\bar{H}_{uq}+\tfrac{1}{2}R$ satisfies

$\bar{H}^{\prime}_{uq}(z)+\left[\bar{H}^{\prime}_{uq}(z)\right]^{H}\geqslant
R,\ \ \forall z\in\mathbb{C}_{\geqslant 1}$

Furthermore, because $A$ is Schur and $R>0$, it follows through a perturbation
argument that for $\rho\in\mathbb{R}_{>0}$ sufficiently small,

$\bar{H}^{\prime}_{uq}((1-\rho)z)+\left[\bar{H}^{\prime}_{uq}((1-\rho)z)\right]^{H}>0,\
\ \forall z\in\mathbb{C}_{\geqslant 1}$

As such, we conclude that $\bar{H}^{\prime}_{uq}(z)$ is strictly positive-real
(SPR). It then follows from the classical Kalman-Yakubovic-Popov Lemma that
there exists a unique $W_{11}=W_{11}^{T}>0$ satisfying (86), with $F_{1}$ as
in (90) and $M>0$ as in (92), and furthermore, that $A_{11}+B_{1}F_{1}$ is
Schur.

Next, consider that (87) can be written in terms of the (known) solution to
$W_{11}$, as

$\displaystyle W_{12}=$
$\displaystyle\left[A_{11}+B_{1}F_{1}\right]^{T}W_{12}A_{22}$
$\displaystyle+A_{11}^{T}W_{11}A_{12}-F_{1}^{T}\left[\tfrac{1}{2}C_{2}-B_{1}^{T}W_{11}A_{12}\right]$

This is a discrete-time Sylvester equation for $W_{12}$, which has a unique
solution because $A_{11}+B_{1}F_{1}$ and $A_{22}$ are Schur. With $W_{12}$
solved, (89) becomes a discrete-time Lyapunov equation for $W_{22}$, which has
a unique solution because $A_{22}$ is Schur.

We conclude that in the case where $(A,C_{q})$ is observable, a unique
$W=W^{T}$ exists which satisfies (35) and renders $A+B_{u}F$ Schur. It remains
only to show that for this observable case, $W>0$. The fact that $W\geqslant
0$ is immediate from (35), because $A$ is Schur and (35) constitutes a
discrete-time Lyapunov equation for $W$ for fixed $F$ and $M$. To show that
$W>0$ it suffices to show that $(A,F)$ is observable. This can be shown by
contradiction; if it were not, there would exist a vector $\eta$ and
eigenvalue $\lambda$ such that $A\eta=\eta\lambda$ and $F\eta=0$. But if this
were the case, it would follow that from a quadratic form on (35) that
$\eta^{H}W\eta\left(1-|\lambda|^{2}\right)=0$. But because $A$ is Schur, it is
known that $|\lambda|<1$ so this implies that $\eta^{H}W\eta=0$. If this is
true, then because $W$ is Hermitian and positive-semidefinite, it must be the
case that $W\eta=0$. But if this is the case, then
$F\eta=-\tfrac{1}{2}M^{-1}C_{q}\eta$. Via the Popov-Belevitch-Hautus test,
this violates the assumption that $(A,C_{q})$ is observable. As such, we
arrive at a contradiction.

###  .3 Proof of Theorem 3

The matrix

$A^{\prime}\triangleq\begin{bmatrix}I&0\\\
B_{u}^{\perp}AB_{u}[C_{d}B_{u}]^{-1}&A_{s}\end{bmatrix}$

is similar to the matrix $A-B_{u}D_{uq}^{-1}C_{q}$. If $\eta$ is an
eigenvector of $A_{s}$, then $[0\ \eta^{H}]^{H}$ is an eigenvector of
$A^{\prime}$. Likewise, the eigenvalues of $A^{\prime}$ are those of $A_{s}$,
together with $\lambda=1$. We conclude that if the eigenvalues of
$A-B_{u}D_{uq}^{-1}C_{q}$ are all in $\mathbb{C}_{<1}\cup\\{1\\}$, then these
same properties hold for $A_{s}$.

Without loss of generality, we assume that the system state space is
partitioned to isolate the subspace controllable from $u$, i.e., we presume

$\displaystyle A=$ $\displaystyle\begin{bmatrix}A_{11}&A_{12}\\\
0&A_{22}\end{bmatrix},$ $\displaystyle B_{u}=$
$\displaystyle\begin{bmatrix}B_{1}\\\ 0\end{bmatrix},$ $\displaystyle C_{q}=$
$\displaystyle\begin{bmatrix}C_{1}^{T}\\\ C_{2}^{T}\end{bmatrix}^{T}$

and we have that

$A-B_{u}D_{uq}^{-1}C_{q}=\begin{bmatrix}\Theta&A_{12}-B_{1}D_{uq}^{-1}C_{2}\\\
0&A_{22}\end{bmatrix}$

where $\Theta\triangleq A_{11}-B_{1}D_{uq}^{-1}C_{1}$. It is known that the
eigenvalues of $A-B_{u}D_{uq}^{-1}C_{q}$ are those of $\Theta$ and $A_{22}$.
Furthermore, the eigenvalues of $A_{22}$ have moduli strictly less than $1$
because $A$ is Schur. So it follows that if $\lambda$ is an eigenvalue of
$A-B_{u}D_{uq}^{-1}C_{q}$ with $|\lambda|\geqslant 1$ then it is also an
eigenvalue of $\Theta$.

Now consider that

$\bar{H}_{uq}^{-1}\sim\left[\begin{array}[]{c|c}\Theta&B_{1}D_{uq}^{-1}\\\
\hline\cr-D_{uq}^{-1}C_{1}&D_{uq}^{-1}\end{array}\right]$

where we note that the above realization is minimal. Consequently,
$\lambda\in\mathbb{C}$ is a pole of $\bar{H}_{uq}^{-1}(z)$ if and only if it
is an eigenvalue of $\Theta$. Because $\bar{H}_{uq}(z)+\bar{H}_{uq}^{H}(z)>0$
for all $z\in\mathbb{C}_{\geqslant 1}\setminus\\{1\\}$, it follows that
$\bar{H}_{uq}(z)$ is nonsingular in this domain, and consequently
$\bar{H}_{uq}^{-1}(z)$ cannot have poles there.

We conclude that all eigenvalues of $A_{s}$ lie in
$\mathbb{C}_{<1}\cup\\{1\\}$. Assume $\lambda=1$ is an eigenvalue of $A_{s}$.
Then there exists a vector $\eta$ such that

$\displaystyle 0=$
$\displaystyle\eta^{H}\left[A_{s}-I\right]=\eta^{H}B_{u}^{\perp}\left[A-I\right]C_{d}^{\perp}[B_{u}^{\perp}C_{d}^{\perp}]^{-1}$

implying that there exists a vector $\nu$ such that

$\eta^{H}B_{u}^{\perp}\left[A-I\right]=\nu^{H}C_{d}.$

Because $A$ is known to be Schur, it follows that $A-I$ is nonsingular, and
consequently $\nu$ must be nonzero. Multiplying from the right by
$[A-I]^{-1}B_{u}$ gives that $\nu^{H}\bar{H}_{ud}(1)=0$. It follows that
$\lambda=1$ is a zero of scalar transfer function $\nu^{H}\bar{H}_{ud}(z)\nu$.
But $\bar{H}_{uq}(z)=\tfrac{1}{T}(z-1)\bar{H}_{ud}(z)$ so we conclude that
$\lambda=1$ must also be a zero of scalar transfer function
$\nu^{H}\bar{H}_{uq}(z)\nu$. If $\nu^{H}\bar{H}_{ud}(z)\nu$ has a zero at
$z=1$, then $\nu^{H}\bar{H}_{uq}(z)\nu$ must have a repeated zero at $z=1$.
But because $\bar{H}_{uq}(z)$ is positive-real, it follows that
$\nu^{H}\bar{H}_{uq}(z)\nu$ is positive-real, precluding the possibility of a
repeated zero at $z=1$. As such, we arrive at a contradiction, and thus
conclude that $\lambda=1$ cannot be an eigenvalue of $A_{s}$.

###  .4 Proof of Theorem 4

We have that for some time $k\in\mathbb{Z}$,

$\displaystyle\operatorname{\mathcal{E}}\\{\bm{p}_{k}\\}=$
$\displaystyle-\mathcal{E}\left\\{\bm{u}^{T}_{k}C_{q}\bm{x}_{k}+\|\bm{u}_{k}\|_{\tilde{R}}^{2}+v_{d}^{T}|\bm{u}_{k}|\right\\}$

where we recall that $\tilde{R}=\tilde{R}^{T}>0$. It is straight-forward to
verify that

$\displaystyle\mathcal{E}\left\\{\|\bm{u}_{k}\|_{\tilde{R}}^{2}\right\\}=$
$\displaystyle\mathcal{E}\left\\{\|\bm{\xi}_{k}\|_{M}^{2}\right\\}+\mathcal{E}\left\\{\|A\bm{x}_{k}+B_{u}\bm{u}_{k}\|_{W}^{2}\right\\}$
$\displaystyle-\mathcal{E}\left\\{\|\bm{x}_{k}\|_{W}^{2}\right\\}-\mathcal{E}\left\\{\bm{u}_{k}^{T}C_{q}\bm{x}_{k}\right\\}$

Because $\bm{x}_{k}$ and $\bm{u}_{k}$ are independent of $\bm{w}_{k}$, and
because $\bm{w}_{k}$ is zero-mean with
$\mathcal{E}\\{\bm{w}_{k}\bm{w}_{k}^{T}\\}=S_{w}$,

$\displaystyle\mathcal{E}\left\\{\|\bm{u}_{k}\|_{\tilde{R}}^{2}\right\\}=$
$\displaystyle\mathcal{E}\left\\{\|\bm{\xi}_{k}\|_{M}^{2}\right\\}-\mathcal{E}\left\\{\bm{u}_{k}^{T}C_{q}\bm{x}_{k}\right\\}-\bar{p}$
$\displaystyle+\mathcal{E}\left\\{\|\bm{x}_{k+1}\|_{W}^{2}\right\\}-\mathcal{E}\left\\{\|\bm{x}_{k}\|_{W}^{2}\right\\}.$

Summing over $k\in\\{0,...,\tau\\}$ and dividing by $\tau+1$ gives the claimed
result, if (54) holds.

###  .5 Proof of Lemma 1

Substitute (70), (71), and (72) into (73), with $J(\cdot,\cdot,\cdot)$ defined
as in (57), to get

$\displaystyle\rho(x_{k+h+1|k})$
$\displaystyle=x_{k+h+1|k}^{T}\left(\sum\limits_{i=0}^{\infty}Q_{i}\right)x_{k+h+1|k}+v_{d}^{T}\left(\sum\limits_{i=0}^{\infty}\left|u_{k+h+1+i|k}\right|\right)$
$\displaystyle\quad+\mu\left(\sum\limits_{i=0}^{\infty}u_{\max}^{T}\max\left\\{0,\left|u_{k+h+1+i|k}\right|-u_{\max}\right\\}\right)$
(93)

where, for $i\in\mathbb{Z}_{\geqslant 0}$,

$\displaystyle Q_{i}\triangleq$
$\displaystyle(B_{u}^{\perp})^{T}(A_{s}^{i})^{T}C_{s}^{T}(F+D_{uq}^{-1}C_{q})^{T}$
$\displaystyle\quad\quad\times
M(F+D_{uq}^{-1}C_{q})C_{s}A_{s}^{i}B_{u}^{\perp}$

By Theorem 3, $A_{s}$ is Schur, and so
$\sum_{i=0}^{\infty}Q_{i}=(B_{u}^{\perp})^{T}ZB_{u}^{\perp}$ where
$Z=Z^{T}\geqslant 0$ is the unique solution to the discrete-time Lyapunov
equation

$Z=A_{s}^{T}ZA_{s}+C_{s}^{T}(F+D_{uq}^{-1}C_{q})^{T}M(F+D_{uq}^{-1}C_{q})C_{s}$

There is no convenient closed-form solution for the second summation in (93),
so we seek to over-bound it by a simple function of $x_{k+h+1|k}$. To do this,
first let $\Pi\in\mathbb{C}^{n\times n}$ be a similarity transformation that
converts $A_{s}$ to its Jordan form $\Lambda$, i.e.,
$A_{s}=\Pi\Lambda\Pi^{-1}$. Then we have that

$\displaystyle\left|u_{k+h+1+i|k}\right|=$
$\displaystyle\sum\limits_{i=0}^{\infty}\left|D_{uq}^{-1}C_{q}C_{s}\Pi\
\Lambda^{i}\ \Pi^{-1}\ B_{u}^{\perp}x_{k+h+1|k}\right|$
$\displaystyle\leqslant$
$\displaystyle\sum\limits_{i=0}^{\infty}\left|D_{uq}^{-1}C_{q}C_{s}\Pi\right|\
\left|\Lambda\right|^{i}\ \left|\Pi^{-1}B_{u}^{\perp}x_{k+h+1|k}\right|$
$\displaystyle=$
$\displaystyle\left|D_{uq}^{-1}C_{q}C_{s}\Pi\right|\left[I-\left|\Lambda\right|\right]^{-1}\left|\Pi^{-1}B_{u}^{\perp}x_{k+h+1|k}\right|$

where we note that the above infinite summation converges because if $\Lambda$
is a Schur Jordan form, then so is $|\Lambda|$.

The last term in (93) can be conservatively over-bounded by recognizing that

$\displaystyle u_{\max}^{T}\max\\{0,|u_{k+h+1+i|k}|-u_{\max}\\}\leqslant$
$\displaystyle\tfrac{1}{4}\|u_{k+h+1+i|k}\|_{2}^{2}$ (94)

But

$\displaystyle\sum\limits_{i=0}^{\infty}\|u_{k+h+1+i|k}\|^{2}=$ $\displaystyle
x_{k+h+1|k}^{T}(B_{u}^{\perp})^{T}XB_{u}^{\perp}x_{k+h+1|k}$

where $X=X^{T}\geqslant 0$ is the solution to Lyapunov equation

$X=A_{s}^{T}XA_{s}+C_{s}^{T}C_{q}^{T}D_{uq}^{-T}D_{uq}^{-1}C_{q}C_{s}.$

We arrive at an over-bound for $\rho(x_{k+h+1|k})$ as in (75), with

$\displaystyle Q_{f}=$ $\displaystyle(B_{u}^{\perp})^{T}(Z+\tfrac{1}{4}\mu
X)B_{u}^{\perp}$ $\displaystyle C_{f}=$
$\displaystyle\operatorname{diag}\left\\{\left[I-|\Lambda|\right]^{-T}\left|D_{uq}^{-1}C_{q}C_{s}\Pi\right|^{T}v_{d}\right\\}\Pi^{T}B_{u}^{\perp}$

To prove inequality (1), consider that because
$x_{k+h+2|k}=C_{s}A_{s}B_{u}^{\perp}x_{k+h+1|k}$ and $B_{u}^{\perp}C_{s}=I$,
we have that

$\displaystyle J_{f}(x_{k+h+2|k})=$
$\displaystyle\left\|B_{u}^{\perp}x_{k+h+1|k}\right\|_{\Xi}^{2}+\|C_{f}x_{k+h+1|k}\|_{1}$

where

$\displaystyle\Xi\triangleq$ $\displaystyle A_{s}^{T}(Z+\tfrac{1}{4}\mu
X)A_{s}$ $\displaystyle=$ $\displaystyle(Z+\tfrac{1}{4}\mu
X)-C_{s}^{T}(F+D_{uq}^{-1}C_{q})^{T}M(F+D_{uq}^{-1}C_{q})C_{s}$
$\displaystyle\quad-\tfrac{1}{4}\mu
C_{s}^{T}C_{q}^{T}D_{uq}^{T}D_{uq}^{-1}C_{q}C_{s}$

and therefore

$\displaystyle\left\|B_{u}^{\perp}x_{k+h+1|k}\right\|_{\Xi}^{2}$
$\displaystyle=\left\|x_{k+h+1|k}\right\|_{Q_{f}}^{2}-\tfrac{1}{4}\mu\left\|u_{k+h+1|k}\right\|^{2}$
$\displaystyle\quad-\left\|u_{k+h+1|k}-Fx_{k+h+1|k}\right\|_{M}^{2}$

Also, consider that

$\displaystyle\|C_{f}x_{k+h+2|k}\|_{1}$
$\displaystyle=v_{d}^{T}\left|D_{uq}^{-1}C_{q}C_{s}\Pi\right|\left[I-|\Lambda|\right]^{-1}\left|\Pi^{-1}B_{u}^{\perp}x_{k+h+2|k}\right|$
$\displaystyle=v_{d}^{T}\left|D_{uq}^{-1}C_{q}C_{s}\Pi\right|\left[I-|\Lambda|\right]^{-1}\left|\Pi^{-1}A_{s}B_{u}^{\perp}x_{k+h+1|k}\right|$
$\displaystyle=v_{d}^{T}\left|D_{uq}^{-1}C_{q}C_{s}\Pi\right|\left[I-|\Lambda|\right]^{-1}\left|\Lambda\Pi^{-1}B_{u}^{\perp}x_{k+h+1|k}\right|$
$\displaystyle\leqslant
v_{d}^{T}\left|D_{uq}^{-1}C_{q}C_{s}\Pi\right|\left[I-|\Lambda|\right]^{-1}|\Lambda|\left|\Pi^{-1}B_{u}^{\perp}x_{k+h+1|k}\right|$
$\displaystyle=\|C_{f}x_{k+h+1|k}\|_{1}$ $\displaystyle\quad-
v_{d}^{T}\left|D_{uq}^{-1}C_{q}C_{s}\Pi\right|\left|\Pi^{-1}B_{u}^{\perp}x_{k+h+1|k}\right|$
$\displaystyle\leqslant\|C_{f}x_{k+h+1|k}\|_{1}-v_{d}^{T}\left|D_{uq}^{-1}C_{q}x_{k+h+1|k}\right|$
$\displaystyle=\|C_{f}x_{k+h+1|k}\|_{1}-v_{d}^{T}\left|u_{k+h+1|k}\right|$

Recalling (94), we conclude (1).

###  .6 Proof of Theorem 5

We begin with two technical lemmas which establish various bounds that will be
useful in the sequel.

###### Lemma 2

Let $K:Y_{k}\mapsto u_{k}$ be as in (64). Then there exist constants
$\\{c_{1},c_{2}\\}\in\mathbb{R}_{>0}$ such that

$\Gamma(x_{k},\mathring{U}_{k|k},\mathring{B}_{k|k})\leqslant
c_{2}\|x_{k}\|_{2}^{2}+c_{1}\|x_{k}\|_{1}$ (95)

######  Proof .1.

For a given $x_{k}$, let the (sub-optimal) sequence $U_{k|k}$ be the sequence
rendering $d_{m|k}=0$ for all $m>k$, i.e.,

$u_{m|k}=\left\\{\begin{array}[]{lll}-\tfrac{1}{T}D_{uq}^{-1}C_{d}Ax_{k}&:&m=k\\\
-D_{uq}^{-1}C_{q}C_{s}A_{s}^{m-k}B_{u}^{\perp}x_{k}&:&m>k\end{array}\right.$

and let corresponding sequence $B_{k|k}$ be the minimum value that satisfies
the force constraint, i.e.,

$b_{m|k}=\max\\{0,|u_{m|k}|-u_{\max}\\}.$

Then the fact that $\\{U_{k|k},B_{k|k}\\}$ are feasible for initial condition
$x_{k}$ implies that that
$\Gamma(x_{k},U_{k|k},B_{k|k})\geqslant\Gamma(x_{k},\mathring{U}_{k|k},\mathring{B}_{k|k})$.
Furthermore, by repeated application of Lemma 1 that

$\displaystyle\Gamma(x_{k},U_{k|k},B_{k|k})$ $\displaystyle\leqslant
J_{f}(x_{k+1|k})+J(x_{k},u_{k|k},b_{k|k})$
$\displaystyle=\|x_{k+1|k}\|_{Q_{f}}^{2}+\|C_{f}x_{k+1|k}\|_{1}+\|u_{k|k}-Fx_{k|k}\|_{M}^{2}$
$\displaystyle\quad+(v_{d}+\mu u_{\max})^{T}|u_{k|k}|$

where we have used the fact that $b_{k|k}\leqslant|u_{k|k}|$. Substitution of
the expression above for $u_{k|k}$, and using the fact that
$C_{s}B_{u}^{\perp}=I$, we have that

$\displaystyle\Gamma(x_{k},U_{k|k},B_{k|k})\leqslant\|\check{A}x_{k}\|_{Q_{f}}^{2}+\|C_{f}\check{A}x_{k}\|_{1}$
$\displaystyle\quad+\|(F+\tfrac{1}{T}D_{uq}^{-1}C_{d})x_{k}\|_{M}^{2}+(v_{d}+\mu
u_{\max})^{T}|\tfrac{1}{T}D_{uq}^{-1}C_{d}x_{k}|$

where $\check{A}\triangleq(I-\tfrac{1}{T}B_{u}D_{uq}^{-1}C_{d})A$. It is then
straight-forward to show by defining $c_{1}$ and $c_{2}$ as

$\displaystyle c_{1}=$
$\displaystyle\|C_{f}\check{A}\|_{1}+\left\|\operatorname{diag}\\{v_{d}+\mu
u_{\max}\\}\tfrac{1}{T}D_{uq}^{-1}C_{d}\right\|_{1}$ $\displaystyle c_{2}=$
$\displaystyle\|Q_{f}^{1/2}\check{A}\|_{2}^{2}+\|M^{1/2}(F+\tfrac{1}{T}D_{uq}^{-1}C_{d})\|_{2}^{2}$

we have that

$\Gamma(x_{k},U_{k|k},B_{k|k})\leqslant
c_{2}\|x_{k}\|_{2}^{2}+c_{1}\|x_{k}\|_{1}.$

Recalling that
$\Gamma(x_{k},U_{k|k},B_{k|k})\geqslant\Gamma(x_{k},\mathring{U}_{k|k},\mathring{B}_{k|k})$
completes the proof.

###### Lemma 3

For two initial states $x_{k}^{i}$, $i\in\\{1,2\\}$, let $U_{k|k}^{i}$ and
$B_{k|k}^{i}$ be feasible sequences for MPC control algorithm (64), i.e.,
trajectories which satisfy

$\displaystyle\Theta(x_{k}^{i},U_{k|k}^{i},B_{k|k}^{i})\leqslant$
$\displaystyle 0,$ $\displaystyle\Psi(x_{k}^{i},U_{k|k}^{i})=$ $\displaystyle
0$ (96)

Let the corresponding receding-horizon trajectories with these input sequences
be $x_{m|k}^{i}$ for $m\geqslant k$. Define perturbations
$\tilde{U}_{k|k}=U_{k|k}^{2}-U_{k|k}^{1}$ and
$\tilde{B}_{k|k}=B_{k|k}^{2}-B_{k|k}^{1}$, and
$\tilde{x}_{m|k}=x_{m|k}^{2}-x_{m|k}^{1}$. Then

$\displaystyle\Gamma(x_{k}^{2},U_{k|k}^{2},B_{k|k}^{2})\leqslant\Gamma(x_{k}^{1},U_{k|k}^{1},B_{k|k}^{1})$
$\displaystyle\quad+2\Gamma^{1/2}(x_{k}^{1},U_{k|k}^{1},B_{k|k}^{1})\chi_{k}^{1/2}+\Gamma(\tilde{x}_{k},\tilde{U}_{k|k},|\tilde{B}_{k|k}|)$
(97)

where

$\chi_{k}\triangleq\sum\limits_{m=k}^{k+h}\|\tilde{u}_{m|k}-F\tilde{x}_{m|k}\|_{M}^{2}+\|\tilde{x}_{k+h+1|k}\|_{Q_{f}}^{2}$
(98)

######  Proof .2.

We have that

$\displaystyle\Gamma(x_{k}^{2},U_{k|k}^{2},B_{k|k}^{2})$
$\displaystyle=\sum\limits_{m=k}^{k+h}\left\|u_{m|k}^{1}-Fx_{m|k}^{1}\right\|_{M}^{2}+\sum\limits_{m=k}^{k+h}\left\|\tilde{u}_{m|k}-F\tilde{x}_{m|k}\right\|_{M}^{2}$
$\displaystyle\quad+2\sum\limits_{m=k}^{k+h}\left(u_{m|k}^{1}-Fx_{m|k}^{1}\right)^{T}M\left(\tilde{u}_{m|k}-F\tilde{x}_{m|k}\right)$
$\displaystyle\quad+\sum\limits_{m=k}^{k+h}\left(v_{d}^{T}\left|u_{m|k}^{1}+\tilde{u}_{m|k}\right|+\mu
u_{\max}^{T}(b_{m|k}^{1}+\tilde{b}_{m|k})\right)$
$\displaystyle\quad+\|x_{k+h+1}^{1}\|_{Q_{f}}^{2}+\|\tilde{x}_{k+h+1}\|_{Q_{f}}^{2}+2\tilde{x}_{k+h+1}^{T}Q_{f}x_{k+h+1}^{1}$
$\displaystyle\quad+\|C_{f}x_{k+h+1}+C_{f}\tilde{x}_{k+h+1}\|_{1}$

Via the Cauchy-Schwartz inequality,

$\displaystyle\sum\limits_{m=k}^{k+h}\left(u_{m|k}^{1}-Fx_{m|k}^{1}\right)^{T}M\left(\tilde{u}_{m|k}-F\tilde{x}_{m|k}\right)$
$\displaystyle+\tilde{x}_{k+h+1}^{T}Q_{f}x_{k+h+1}^{1}$
$\displaystyle\leqslant\left(\sum\limits_{m=k}^{k+h}\left\|u_{m|k}^{1}-Fx_{m|k}^{1}\right\|_{M}^{2}+\|x_{k+h+1}^{1}\|_{Q_{f}}^{2}\right)^{1/2}\chi_{k}$
$\displaystyle\leqslant\Gamma^{1/2}(x_{k}^{1},U_{k|k}^{1},B_{k|k}^{1})\chi_{k}$

and via the triangle inequality,

$\displaystyle\sum\limits_{m=k}^{k+h}v_{d}^{T}\left|u_{m|k}^{1}+\tilde{u}_{m|k}\right|+\left\|C_{f}x_{k+h+1}^{1}+C_{f}\tilde{x}_{k+h+1}\right\|_{1}$
$\displaystyle\leqslant\sum\limits_{m=k}^{k+h}v_{d}^{T}\left|u_{m|k}^{1}\right|+\sum\limits_{m=k}^{k+h}v_{d}^{T}\left|\tilde{u}_{m|k}\right|$
$\displaystyle\quad+\left\|C_{f}x_{k+h+1}^{1}\right\|_{1}+\left\|C_{f}\tilde{x}_{k+h+1}\right\|_{1}$

Also, note that because $b_{m|k}^{1}$ and $b_{m|k}^{2}$ are both nonnegative,

$u_{\max}^{T}(b_{m|k}^{1}+\tilde{b}_{m|k})\leqslant
u_{\max}^{T}b_{m|k}^{1}+u_{\max}^{T}|\tilde{b}_{m|k}|.$

Substituting these inequalities gives (97).

Let the present time be $k$, and let $\mathring{U}_{k|k}$ and
$\mathring{B}_{k|k}$ be the trajectories found by the MPC algorithm (64),
given present state $x_{k}$. For $m\geqslant k$, let the optimized receding-
horizon state trajectory be denoted $\mathring{x}_{m|k}$, with
$\mathring{x}_{k|k}=x_{k}$. Let $\mathring{d}_{m|k}$ be as in (68), evaluated
with the optimal trajectories, i.e.,
$\mathring{d}_{m|k}=C_{d}\mathring{x}_{m|k}$, where we note that
$\mathring{d}_{k+h+1|k}=0$. For $m>k+h$, we denote the post-horizon input
$\mathring{u}_{m|k}$ as the input that maintains $\mathring{d}_{m|k}=0$, i.e.,

$\displaystyle\mathring{x}_{m|k}=$ $\displaystyle
C_{s}A_{s}^{m-(k+h+1)}B_{u}^{\perp}\mathring{x}_{k+h+1|k}$
$\displaystyle\mathring{u}_{m|k}=$ $\displaystyle-
D_{uq}^{-1}C_{q}\mathring{x}_{m|k}$

Now, we advance the present time to $k+1$. With acquisition of data $y_{k}$,
the data $w_{k}$ and $x_{k+1}$ become known. We consider the receding-horizon
control trajectory $U_{k+1|k+1}$ that enforces

$d_{m|k+1}=\mathring{d}_{m|k},\ \ \forall m\in\\{k+2,...,k+h+2\\}$ (99)

We formulate the corresponding receding-horizon slack variable trajectory
$B_{k+1|k+1}$ as

$b_{m|k+1}=\max\left\\{0,|u_{m|k+1}|-u_{\max}\right\\}$

Note that, so formulated, this trajectory pair $\\{U_{k+1|k+1},B_{k+1|k+1}\\}$
satisfies the constraints in the MPC algorithm at the $(k+1)^{th}$ iteration,
i.e., conditions

$\displaystyle\Theta(x_{k+1},U_{k+1|k+1},B_{k+1|k+1})\leqslant$ $\displaystyle
0,$ $\displaystyle\Psi(x_{k+1},U_{k+1|k+1})=$ $\displaystyle 0$

and is therefore in the feasibility domain for the MPC optimization algorithm
at time $k+1$.

We now claim that the sequence $U_{k+1|k+1}$ resulting in (99), and the
associated receding-horizon sequence $\\{x_{m|k+1}:m\geqslant k+1\\}$, can be
expressed as

$\displaystyle u_{m|k+1}=$ $\displaystyle\mathring{u}_{m|k}+\tilde{u}_{m|k+1}$
(100) $\displaystyle x_{m|k+1}=$
$\displaystyle\mathring{x}_{m|k}+\tilde{x}_{m|k+1}$ (101)

where

$\displaystyle\tilde{u}_{m|k+1}$
$\displaystyle=-\tfrac{1}{T}D_{uq}^{-1}C_{d}A\tilde{x}_{m|k+1}$ (102)
$\displaystyle\tilde{x}_{m|k+1}$
$\displaystyle=\left\\{\begin{array}[]{lll}B_{w}w_{k}&:&m=k+1\\\
\tilde{A}^{m-k-1}B_{w}w_{k}&:&m>k+1\end{array}\right.$ (105)

with $\tilde{A}\triangleq
C_{d}^{\perp}[B_{u}^{\perp}C_{d}^{\perp}]^{-1}B_{u}^{\perp}A$. To show this,
consider the inverse system dynamic model from Section 2.4. Let

$\displaystyle\zeta_{m|k+1}\triangleq$
$\displaystyle\mathcal{E}\\{\bm{\zeta}_{m}|\bm{Y}_{k+1}=Y_{k+1},\
\bm{d}_{m}=d_{m|k+1},\forall m\geqslant k+1\\}.$

Then we have that for $m\geqslant k+1$,

$\displaystyle\zeta_{m+1|k+1}=$ $\displaystyle
A_{s}\zeta_{m|k+1}+B_{ds}d_{m|k+1}$

with initial condition

$\zeta_{k+1|k+1}=A_{s}\zeta_{k}+B_{ws}w_{k}+B_{ds}d_{k}$

But we also have that for $m\geqslant k+1$,

$\displaystyle\mathring{\zeta}_{m+1|k}=$ $\displaystyle
A_{s}\mathring{\zeta}_{m|k}+B_{ds}\mathring{d}_{m|k}$

with initial condition

$\mathring{\zeta}_{k+1|k}=A_{s}\zeta_{k}+B_{ds}\mathring{d}_{k|k}$

So, letting
$\tilde{\zeta}_{m|k+1}\triangleq\zeta_{m|k+1}-\mathring{\zeta}_{m|k}$, we have
that

$\tilde{\zeta}_{m+1|k+1}=A_{s}\tilde{\zeta}_{m|k+1}+B_{ds}\left(d_{m|k+1}-\mathring{d}_{m|k}\right)$

with initial condition

$\tilde{\zeta}_{k+1|k+1}=B_{ws}w_{k}+B_{ds}\left(d_{k}-\mathring{d}_{k|k}\right)$

But $d_{k}-\mathring{d}_{k|k}=D_{wd}w_{k}$, so using (45), the above initial
condition simplifies to

$\tilde{\zeta}_{k+1|k+1}=B_{u}^{\perp}B_{w}w_{k}$

Also, we have that $d_{k+1|k+1}-\mathring{d}_{k+1|k}=C_{d}B_{w}w_{k}$ while
for $m>k+1$, (99) implies that $d_{m|k+1}-\mathring{d}_{m|k}=0.$ Consequently,
we have that for $m>k+1$,

$\displaystyle\tilde{\zeta}_{m|k+1}=$ $\displaystyle
A_{s}^{m-k-2}\left[A_{s}B_{u}^{\perp}+B_{ds}C_{d}\right]B_{w}w_{k}$
$\displaystyle=$ $\displaystyle A_{s}^{m-k-2}B_{u}^{\perp}AB_{w}w_{k}$

where we have used (44) and (45) in the second line. For $m\geqslant k+1$, it
follows from (43) and the above, that

$\displaystyle\tilde{x}_{m|k+1}=$
$\displaystyle\left\\{\begin{array}[]{lll}C_{s}B_{u}^{\perp}B_{w}w_{k}+D_{ds}C_{d}B_{w}w_{k}&:&m=k+1\\\
C_{s}A_{s}^{m-k-2}B_{u}^{\perp}AB_{w}w_{k}&:&m>k+1\end{array}\right.$

which simplifies to (105). Using (38), this implies that

$\displaystyle\tilde{u}_{m|k+1}=-D_{uq}^{-1}C_{q}\tilde{x}_{m|k+1}+\tfrac{1}{T}D_{uq}^{-1}\left(\tilde{d}_{m+1|k+1}-\tilde{d}_{m|k+1}\right).$

Substitution of (105) and use of (44)-(47) give the above as equivalent to
(102).

Define sequences
$\mathring{U}_{k+1|k}=\\{\mathring{u}_{k+1|k},...,\mathring{u}_{k+h+1|k}\\}$
and
$\mathring{B}_{k+1|k}=\\{\mathring{b}_{k+1|k},...,\mathring{b}_{k+h+1|k}\\}$.
Then we have that the pair $\\{\mathring{U}_{k+1|k},\mathring{B}_{k+1|k}\\}$
are in the feasibility domain of MPC optimization algorithm (64) at time $k+1$
given initial condition $\mathring{x}_{k+1|k}$. From Lemma 3, we infer that

$\displaystyle\Gamma(x_{k+1},U_{k+1|k+1},B_{k+1|k+1})$
$\displaystyle\leqslant\Gamma(\mathring{x}_{k+1|k},\mathring{U}_{k+1|k},\mathring{B}_{k+1|k})$
$\displaystyle\quad+2\Gamma^{1/2}(\mathring{x}_{k+1|k},\mathring{U}_{k+1|k},\mathring{B}_{k+1|k})\chi_{k+1}$
$\displaystyle\quad+\Gamma(\tilde{x}_{k+1},\tilde{U}_{k+1|k+1},|\tilde{U}_{k+1|k+1}|)$

where

$\chi_{k+1}\triangleq\sum\limits_{m=k+1}^{k+h+1}\|\tilde{u}_{m|k+1}-F\tilde{x}_{m|k+1}\|_{M}^{2}+\|\tilde{x}_{k+h+2|k+1}\|_{Q_{f}}^{2}$
(106)

From Lemma 1 we have inequality (1), and therefore

$\Gamma(x_{k},\mathring{U}_{k|k},\mathring{B}_{k|k})\geqslant
J(x_{k},u_{k},b_{k})+\Gamma(\mathring{x}_{k+1|k},\mathring{U}_{k+1|k},\mathring{B}_{k+1|k}).$

So it follows that

$\Gamma(x_{k+1},U_{k+1|k+1},B_{k+1|k+1})\leqslant\Gamma(x_{k},\mathring{U}_{k|k},\mathring{B}_{k|k})\\\
-J(x_{k},u_{k},b_{k})+2\Gamma^{1/2}(\mathring{x}_{k+1|k},\mathring{U}_{k+1|k},\mathring{B}_{k+1|k})\chi_{k+1}\\\
+\Gamma(\tilde{x}_{k},\tilde{U}_{k|k},|\tilde{U}_{k|k}|)$

Now, using Lemma 2, we note that

$\displaystyle\Gamma^{1/2}(\mathring{x}_{k+1|k},\mathring{U}_{k+1|k},\mathring{B}_{k+1|k})$
$\displaystyle\leqslant\left[\Gamma(\mathring{x}_{k|k},\mathring{U}_{k|k},\mathring{B}_{k|k})-J(x_{k},u_{k},b_{k})\right]^{1/2}$
$\displaystyle\leqslant\Gamma^{1/2}(\mathring{x}_{k|k},\mathring{U}_{k|k},\mathring{B}_{k|k})$
$\displaystyle\leqslant\left[c_{2}\|x_{k}\|_{2}^{2}+c_{1}\|x_{k}\|_{1}\right]^{1/2}$
$\displaystyle\leqslant\left[c_{2}\|x_{k}\|_{2}^{2}+c_{1}\sqrt{n}\|x_{k}\|_{2}\right]^{1/2}$
$\displaystyle\leqslant\sqrt{c_{2}}\left[\|x_{k}\|_{2}+\frac{c_{1}\sqrt{n}}{2c_{2}}\right]$

So we conclude that

$\Gamma(x_{k+1},U_{k+1|k+1},B_{k+1|k+1})\leqslant\Gamma(x_{k},\mathring{U}_{k|k},\mathring{B}_{k|k})\\\
-J(x_{k},u_{k},b_{k})+2\sqrt{c_{2}}\left[\|x_{k}\|_{2}+\frac{c_{1}\sqrt{n}}{2c_{2}}\right]\chi_{k+1}\\\
+\Gamma(\tilde{x}_{k},\tilde{U}_{k|k},|\tilde{U}_{k|k}|)$

Now, we recognize that due to (102) and (105), together with (66) and (106),
it follows that there exist constants $a_{1}$, $a_{2}$, and $a_{3}$ such that

$\displaystyle\chi_{k+1}^{1/2}\leqslant a_{1}\|w_{k}\|_{2}\leqslant
a_{1}\|w_{k}\|_{1}$
$\displaystyle\Gamma(\tilde{x}_{k},\tilde{U}_{k|k},|\tilde{U}|_{k|k})\leqslant
a_{2}\|w_{k}\|_{2}^{2}+a_{3}\|w_{k}\|_{1}$

As such, we conclude that

$\displaystyle J(x_{k},u_{k},b_{k})$
$\displaystyle\leqslant\Gamma(x_{k},\mathring{U}_{k|k},\mathring{B}_{k|k})-\Gamma(x_{k+1},U_{k+1|k+1},B_{k+1|k+1})$
$\displaystyle\quad+a_{2}\|w_{k}\|_{2}^{2}+a_{4}\|w_{k}\|_{1}+a_{5}\|x_{k}\|_{2}\|w_{k}\|_{2}$

where $a_{4}\triangleq a_{3}+a_{1}c_{1}\sqrt{n/c_{2}}$ and $a_{5}\triangleq
2\sqrt{c_{2}}a_{1}$.

Now consider the implementation of MPC optimization algorithm (64) at time
$k+1$, resulting in optimal sequences $\mathring{U}_{k+1|k+1}$ and
$\mathring{B}_{k+1|k+1}$. Because sequences $U_{k+1|k+1}$ and $B_{k+1|k+1}$
derived above are feasible, it follows that

$\Gamma(x_{k+1},\mathring{U}_{k+1|k+1},\mathring{B}_{k+1|k+1})\\\
\leqslant\Gamma(x_{k+1},U_{k+1|k+1},B_{k+1|k+1})$

Define $r_{\tau}$ for all $\tau\in\mathbb{Z}_{\geqslant 0}$ as

$r_{\tau}\triangleq\sum\limits_{k=0}^{\tau}J(x_{k},u_{k},b_{k})$

Then the above inequalities imply that

$\displaystyle r_{\tau}\leqslant$
$\displaystyle\Gamma(x_{0},\mathring{U}_{0|0},\mathring{B}_{0|0})+a_{2}\|w\|_{2\tau}^{2}+a_{4}\|w\|_{1\tau}$
$\displaystyle+a_{5}\sum\limits_{k=0}^{\tau}\|x_{k}\|_{2}\|w_{k}\|_{2}$

Using Lemma 2 and the Cauchy-Schwartz inequality,

$\displaystyle r_{\tau}\leqslant$ $\displaystyle
c_{2}\|x_{0}\|_{2}^{2}+c_{1}\|x_{0}\|_{1}+a_{2}\|w\|_{2\tau}^{2}+a_{4}\|w\|_{1\tau}$
$\displaystyle+a_{5}\|x\|_{2\tau}\|w\|_{2\tau}$

Let $\delta_{k}\triangleq M^{1/2}(u_{k}-Fx_{k})$. Then it follows that
$r_{\tau}\geqslant\|\delta\|_{2\tau}^{2}$, and thus we have that

$\displaystyle\|\delta\|_{2\tau}^{2}\leqslant$ $\displaystyle
c_{2}\|x_{0}\|_{2}^{2}+c_{1}\|x_{0}\|_{1}+a_{2}\|w\|_{2\tau}^{2}+a_{4}\|w\|_{1\tau}$
$\displaystyle+a_{5}\|x\|_{2\tau}\|w\|_{2\tau}$

Meanwhile, consider that

$x_{k+1}=\left[A+B_{u}F\right]x_{k}+B_{u}M^{-1/2}\delta_{k}+B_{w}w_{k}$

where we recall from Theorem 2 that $[A+B_{u}F]$ is Schur. It follows that
there exist $\\{\vartheta,\varphi\\}\subset\mathbb{R}_{>0}$ such that

$\|x\|_{2\tau}^{2}\leqslant\vartheta\|\delta\|_{2\tau}^{2}+\varphi\|w\|_{2\tau}^{2}$

for all $\tau\in\mathbb{R}_{>0}$. So we conclude that

$\displaystyle 0$ $\displaystyle\geqslant\|x\|_{2\tau}^{2}-\vartheta
a_{5}\|w\|_{2\tau}\|x\|_{2\tau}$ $\displaystyle-\left(\vartheta
c_{2}\|x_{0}\|_{2}^{2}+\vartheta c_{1}\|x_{0}\|_{1}+(\vartheta
a_{2}+\varphi)\|w\|_{2\tau}^{2}+\vartheta a_{4}\|w\|_{1\tau}\right)$

The right-hand side is convex in $\|x\|_{2\tau}$, with real roots. So

$\|x\|_{2\tau}\leqslant\tfrac{1}{2}\vartheta a_{5}\|w\|_{2\tau}\\\
+\sqrt{\begin{array}[]{l}\vartheta c_{2}\|x_{0}\|_{2}^{2}+\vartheta
c_{1}\|x_{0}\|_{1}+\vartheta a_{4}\|w\|_{1\tau}\\\ +(\vartheta
a_{2}+\varphi+\tfrac{1}{4}\vartheta^{2}a_{5}^{2})\|w\|_{2\tau}^{2}\end{array}}$

Inequality (78) is obtained by squaring the above expression and linearizing
the (concave) square-root term about the value $\tfrac{1}{2}\vartheta
a_{5}\|w\|_{2\tau}$.

## References

  * [1] M. E. McCormick, _Ocean wave energy conversion_. Courier Corporation, 2013.
  * [2] F. d. O. Antonio, “Wave energy utilization: A review of the technologies,” _Renewable and sustainable energy reviews_ , vol. 14, no. 3, pp. 899–918, 2010.
  * [3] D. Clemente, P. Rosa-Santos, and F. Taveira-Pinto, “On the potential synergies and applications of wave energy converters: A review,” _Renewable and Sustainable Energy Reviews_ , vol. 135, p. 110162, 2021.
  * [4] B. Guo and J. V. Ringwood, “A review of wave energy technology from a research and commercial perspective,” _IET Renewable Power Generation_ , vol. 15, no. 14, pp. 3065–3090, 2021.
  * [5] J. Falnes, “A review of wave-energy extraction,” _Marine structures_ , vol. 20, no. 4, pp. 185–201, 2007.
  * [6] S. H. Salter, J. Taylor, and N. Caldwell, “Power conversion mechanisms for wave energy,” _Proceedings of the Institution of Mechanical Engineers, Part M: Journal of Engineering for the Maritime Environment_ , vol. 216, no. 1, pp. 1–27, 2002.
  * [7] M. A. Jusoh, M. Z. Ibrahim, M. Z. Daud, A. Albani, and Z. Mohd Yusop, “Hydraulic power take-off concepts for wave energy conversion system: A review,” _Energies_ , vol. 12, no. 23, p. 4510, 2019.
  * [8] R. Ahamed, K. McKee, and I. Howard, “Advancements of wave energy converters based on power take off (pto) systems: A review,” _Ocean Engineering_ , vol. 204, p. 107248, 2020.
  * [9] ——, “A review of the linear generator type of wave energy converters? power take-off systems,” _Sustainability_ , vol. 14, no. 16, p. 9936, 2022\.
  * [10] D. B. Murray, J. G. Hayes, D. L. O’Sullivan, and M. G. Egan, “Supercapacitor testing for power smoothing in a variable speed offshore wave energy converter,” _IEEE Journal of Oceanic Engineering_ , vol. 37, no. 2, pp. 301–308, 2012.
  * [11] E. Tedeschi and M. Santos-Mugica, “Modeling and control of a wave energy farm including energy storage for power quality enhancement: The bimep case study,” _IEEE Transactions on Power Systems_ , vol. 29, no. 3, pp. 1489–1497, 2013.
  * [12] J. Hals, J. Falnes, and T. Moan, “A comparison of selected strategies for adaptive control of wave energy converters,” _Journal of Offshore Mechanics and Arctic Engineering_ , vol. 133, no. 3, 2011.
  * [13] J. Scruggs, S. Lattanzio, A. Taflanidis, and I. Cassidy, “Optimal causal control of a wave energy converter in a random sea,” _Applied Ocean Research_ , vol. 42, pp. 1–15, 2013.
  * [14] J. V. Ringwood, G. Bacelli, and F. Fusco, “Energy-maximizing control of wave-energy converters: The development of control system technology to optimize their operation,” _IEEE control systems magazine_ , vol. 34, no. 5, pp. 30–55, 2014.
  * [15] R. G. Coe, G. Bacelli, D. G. Wilson, O. Abdelkhalik, U. A. Korde, and R. D. Robinett III, “A comparison of control strategies for wave energy converters,” _International journal of marine energy_ , vol. 20, pp. 45–63, 2017.
  * [16] J. Falnes, “Radiation impedance matrix and optimum power absorption for interacting oscillators in surface waves,” _Applied ocean research_ , vol. 2, no. 2, pp. 75–80, 1980.
  * [17] D. Evans, “Power from water waves,” _Annual review of Fluid mechanics_ , vol. 13, no. 1, pp. 157–187, 1981.
  * [18] P. Nebel, “Maximizing the efficiency of wave-energy plant using complex-conjugate control,” _Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering_ , vol. 206, no. 4, pp. 225–236, 1992.
  * [19] S. Zhan and G. Li, “Linear optimal noncausal control of wave energy converters,” _IEEE Transactions on Control Systems Technology_ , vol. 27, no. 4, pp. 1526–1536, 2018.
  * [20] D. Evans, “Maximum wave-power absorption under motion constraints,” _Applied Ocean Research_ , vol. 3, no. 4, pp. 200–203, 1981.
  * [21] J. Hals, J. Falnes, and T. Moan, “Constrained optimal control of a heaving buoy wave-energy converter,” _Journal of Offshore Mechanics and Arctic Engineering_ , vol. 133, no. 1, 2011.
  * [22] G. Bacelli and J. Ringwood, “Constrained control of arrays of wave energy devices,” _International Journal of Marine Energy_ , vol. 3, pp. e53–e69, 2013.
  * [23] M. Penalba, G. Giorgi, and J. V. Ringwood, “Mathematical modelling of wave energy converters: A review of nonlinear approaches,” _Renewable and Sustainable Energy Reviews_ , vol. 78, pp. 1188–1207, 2017.
  * [24] J. Cretel, A. Lewis, G. Lightbody, and G. Thomas, “An application of model predictive control to a wave energy point absorber,” _IFAC Proceedings Volumes_ , vol. 43, no. 1, pp. 267–272, 2010.
  * [25] M. Richter, M. E. Magana, O. Sawodny, and T. K. Brekken, “Nonlinear model predictive control of a point absorber wave energy converter,” _IEEE Transactions on Sustainable Energy_ , vol. 4, no. 1, pp. 118–126, 2012.
  * [26] G. Li and M. R. Belmont, “Model predictive control of sea wave energy converters–part i: A convex approach for the case of a single device,” _Renewable Energy_ , vol. 69, pp. 453–463, 2014.
  * [27] N. Tom and R. W. Yeung, “Experimental confirmation of nonlinear-model-predictive control applied offline to a permanent magnet linear generator for ocean-wave energy conversion,” _IEEE Journal of Oceanic Engineering_ , vol. 41, no. 2, pp. 281–295, 2015.
  * [28] N. Faedo, S. Olaya, and J. V. Ringwood, “Optimal control, mpc and mpc-like algorithms for wave energy systems: An overview,” _IFAC Journal of Systems and Control_ , vol. 1, pp. 37–56, 2017.
  * [29] Q. Zhong and R. W. Yeung, “An efficient convex formulation for model-predictive control on wave-energy converters,” _Journal of Offshore Mechanics and Arctic Engineering_ , vol. 140, no. 3, 2018.
  * [30] F. Allgöwer and A. Zheng, _Nonlinear model predictive control_. Birkhäuser, 2012, vol. 26.
  * [31] G. Li, “Nonlinear model predictive control of a wave energy converter based on differential flatness parameterisation,” _International Journal of Control_ , vol. 90, no. 1, pp. 68–77, 2017.
  * [32] J. Na, B. Wang, G. Li, S. Zhan, and W. He, “Nonlinear constrained optimal control of wave energy converters with adaptive dynamic programming,” _IEEE Transactions on Industrial Electronics_ , vol. 66, no. 10, pp. 7904–7915, 2018.
  * [33] A. S. Haider, T. K. Brekken, and A. McCall, “Real-time nonlinear model predictive controller for multiple degrees of freedom wave energy converters with non-ideal power take-off,” _Journal of Marine Science and Engineering_ , vol. 9, no. 8, p. 890, 2021.
  * [34] M. Previsic, A. Karthikeyan, and D. Lyzenga, “In-ocean validation of a deterministic sea wave prediction (dswp) system leveraging x-band radar to enable optimal control in wave energy conversion systems,” _Applied Ocean Research_ , vol. 114, p. 102784, 2021.
  * [35] O. Abdelkhalik, S. Zou, R. Robinett, G. Bacelli, and D. Wilson, “Estimation of excitation forces for wave energy converters control using pressure measurements,” _International Journal of Control_ , vol. 90, no. 8, pp. 1793–1805, 2017.
  * [36] F. Fusco and J. V. Ringwood, “Short-term wave forecasting for real-time control of wave energy converters,” _IEEE Transactions on sustainable energy_ , vol. 1, no. 2, pp. 99–106, 2010.
  * [37] M. P. Schoen, J. Hals, and T. Moan, “Wave prediction and robust control of heaving wave energy devices for irregular waves,” _IEEE Transactions on energy conversion_ , vol. 26, no. 2, pp. 627–638, 2011.
  * [38] G. Li, G. Weiss, M. Mueller, S. Townley, and M. R. Belmont, “Wave energy converter control by wave prediction and dynamic programming,” _Renewable Energy_ , vol. 48, pp. 392–403, 2012.
  * [39] U. A. Korde, “Near-optimal control of a wave energy device in irregular waves with deterministic-model driven incident wave prediction,” _Applied Ocean Research_ , vol. 53, pp. 31–45, 2015.
  * [40] A. F. Davis and B. C. Fabien, “Wave excitation force estimation of wave energy floats using extended kalman filters,” _Ocean Engineering_ , vol. 198, p. 106970, 2020.
  * [41] S. Mavrakos and P. McIver, “Comparison of methods for computing hydrodynamic characteristics of arrays of wave power devices,” _Applied Ocean Research_ , vol. 19, no. 5-6, pp. 283–291, 1997.
  * [42] M. Folley, A. Babarit, B. Child, D. Forehand, L. O?Boyle, K. Silverthorne, J. Spinneken, V. Stratigaki, and P. Troch, “A review of numerical modelling of wave energy converter arrays,” in _International Conference on Offshore Mechanics and Arctic Engineering_ , vol. 44946. American Society of Mechanical Engineers, 2012, pp. 535–545.
  * [43] Y. Lao, J. T. Scruggs, A. Karthikeyan, and M. Previsic, “Discrete-time causal control of a wave energy converter with finite stroke in stochastic waves,” _IEEE Transactions on Control Systems Technology_ , vol. 30, no. 3, pp. 1198–1214, 2022.
  * [44] A. H. Sayed and T. Kailath, “A survey of spectral factorization methods,” _Numerical Linear Algebra with Applications_ , vol. 8, no. 6-7, pp. 467–496, 2001.
  * [45] B. Anderson and J. Moore, _Optimal Filtering_ , ser. Dover Books on Electrical Engineering. Dover Publications, 2012. [Online]. Available: https://books.google.com/books?id=iYMqLQp49UMC
  * [46] J. Hoagg, S. Lacy, R. Erwin, and D. Bernstein, “First-order-hold sampling of positive real systems and subspace identification of positive real models,” in _Proceedings of the 2004 American Control Conference_ , vol. 1, 2004, pp. 861–866 vol.1.
  * [47] Y. Lao and J. T. Scruggs, “A modified technique for spectral factorization of infinite-dimensional systems using subspace techniques,” in _2019 IEEE 58th Conference on Decision and Control (CDC)_ , 2019, pp. 5412–5419.
  * [48] O. Faltinsen, _Sea loads on ships and offshore structures_. Cambridge university press, 1993, vol. 1.
  * [49] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex programming, version 2.1,” http://cvxr.com/cvx, Mar. 2014.

[]Connor Ligeikis (S’21) is a Ph.D. Candidate in the Department of Civil and
Environmental Engineering at the University of Michigan. He received his
B.S.E. and M.S. degrees in Civil Engineering from the University of
Connecticut in 2017 and 2019, respectively, and his M.S. in Electrical and
Computer Engineering from the University of Michigan in 2021. His research
interests include vibration, control, mechatronics, and cyber-physical
systems. He is a National Science Foundation Graduate Research Fellow.

[]Jeff Scruggs (M’11) is an Associate Professor in the Departments of Civil
and Environmental Engineering, as well as of Electrical and Computer
Engineering (by courtesy), at the University of Michigan, which he joined in
2011. He received his B.S. and M.S. degrees in Electrical Engineering from
Virginia Tech in 1997 and 1999, respectively, and his Ph.D. in Applied
Mechanics from the Caltech in 2004. Prior to joining the University of
Michigan, he held postdoctoral positions at Caltech and the University of
California, San Diego, and was on the faculty at Duke University from 2007-11.
Scruggs’s current research is in the areas of mechanics, vibration, energy,
and control.