# Nonlinear coupling of reversed shear Alfvén eigenmode and toroidal Alfvén
eigenmode during current ramp

Shizhao Wei1, Yahui Wang1, Peiwan Shi2,3, Wei Chen2, Ningfei Chen1 and Zhiyong
Qiu1111E-mail<EMAIL_ADDRESS>1Institute for Fusion Theory and Simulation and
Department of Physics, Zhejiang University, Hangzhou 310027, P.R.C
2Southwestern Institute of Physics - P.O. Box 432 Chengdu 610041, P.R.C.
3Key Laboratory of Materials Modification by Laser, Ion, and Electron Beams
(Ministry of Education), School of Physics, Dalian University of Technology,
Dalian 116024, P.R.C.

###### Abstract

Two novel nonlinear mode coupling processes for reversed shear Alfvén
eigenmode (RSAE) nonlinear saturation are proposed and investigated. In the
first process, RSAE nonlinearly couples to a co-propagating toroidal Alfvén
eigenmode (TAE) with the same toroidal and poloidal mode numbers, and
generates a geodesic acoustic mode (GAM). In the second process, RSAE couples
to a counter-propagating TAE and generates an ion acoustic wave quasi-mode
(IAW). The condition for the two processes to occur is favored during current
ramp. Both processes contribute to effectively saturate the Alfvénic
instabilities, as well as nonlinearly transfer of energy from energetic fusion
alpha particles to fuel ions in burning plasmas.

In next generation magnetically confined fusion devices such as ITER Tomabechi
et al. (1991), energetic particles (EPs), e.g., fusion alpha particles, are
expected to play important roles as they contribute significantly to the total
power density and drive symmetry breaking collective modes including shear
Alfvén wave (SAW) instabilities Fasoli et al. (2007); Chen and Zonca (2016);
Chen and Wang (2020). SAW can lead to EP anomalous transport loss, degradation
of burning plasma performance and possible damage of plasma facing components
Ding et al. (2015). So for understanding of EP confinement and thus fusion
plasma performance, the in-depth research of SAW instability dynamics
including nonlinear evolution is needed. Due to plasma nonuniformities and
equilibrium magnetic field geometry, SAW instabilities can be excited as EP
modes (EPMs) in the continuum Chen (1994), or various discrete Alfvén
eigenmodes (AEs) inside forbidden gaps of SAW continuum such as toroidal AE
(TAE) Cheng et al. (1985); Chen (1988); Fu and Van Dam (1989), reversed shear
AE (RSAE) Berk et al. (2001); Zonca et al. (2002); Kramer et al. (2004), beta-
induced AE (BAE) Heidbrink et al. (1993), etc. The nonlinear mode coupling of
SAW instabilities, providing an effective channel for the mode saturation, has
been observed in experiments Buratti et al. (2005), and investigated
analytically Hahm and Chen (1995); Zonca et al. (1995); Todo et al. (2010);
Biancalani et al. (2010); Chen and Zonca (2012) as well as in large scale
simulations Todo et al. (2010); Zhang and Lin (2013); Biancalani et al.
(2018). In this work, two possible channels for RSAE saturation as the result
of the nonlinear coupling to TAE during current ramp are proposed and
investigated, which are of relevance for burning plasmas in future reactors
due to their typically reversed shear advanced scenarios where RSAEs can be
strongly driven unstable by core localized fusion alpha particles Tomabechi et
al. (1991); Gong et al. (2019); Huang et al. (2020); Wang et al. (2018); Ren
et al. (2020).

TAE is one of the most well-known SAW eigenmodes and widely exists in present-
day and future magnetically confined plasmas. It can be resonantly excited by
EPs in the SAW continuum gap induced by torodicity Cheng et al. (1985); Fu and
Van Dam (1989), as the result of the coupling of neighbouring poloidal
harmonics of SAW continuum, and is typically localized near the center of two
mode rational surfaces, where
$q\left(r\right)\simeq\left(2m+1\right)/\left(2n\right)$, with the parallel
wavenumbers of the two dominant poloidal harmonics both being
$\left|k_{\parallel}\right|\simeq 1/\left(2qR_{0}\right)$. Here, $m/n$ are the
poloidal/toroidal mode number and $q$ is the safety factor of the torus.
Furthermore, TAE is considered as an important player in transporting EPs due
to its low excitation threshold and suitable resonance condition with fusion
alpha particles with $v_{\alpha}\gtrsim v_{A}$ in ITER Fasoli et al. (2007).
Here, $v_{\alpha}$ is the birth velocity of the fusion alpha particle, and
$v_{A}$ is the Alfvén speed.

RSAE, also called Alfvén cascade in literatures, is excited near the minima of
the safety factor ($q_{min}$), and is often composed by one dominant poloidal
harmonic Berk et al. (2001); Zonca et al. (2002). With the parallel wave
number $k_{\parallel}\simeq\left|n-m/q_{min}\right|/R$ and thus its frequency
$\omega\simeq k_{\parallel}v_{A}$ determined by the value of $q_{min}$, RSAE
related physics is thus, sensitive to the $q$-profile. One thus expects that,
with the change of $q$-profile during, e.g., current ramp, the RSAE frequency
can sweep from BAE (as $\left|nq_{min}-m\right|\simeq 0$) to TAE (as
$\left|nq_{min}-m\right|\simeq 1/2$) frequency range, which has been shown in
many experiments Van Zeeland et al. (2007) as well as simulations Wang et al.
(2020). Furthermore, it is shown in Ref. Van Zeeland et al. (2007) that, as
the core-localized RSAE frequency sweeps up, RSAE can temporally couple to TAE
and generate a global mode, creating an effective channel for particle global
transport from Tokamak core to edge. This channel was also predicted and
investigated theoretically in Ref. Zonca et al. (2002). The coupling of RSAE
and TAE with different mode numbers and generating a low frequency mode during
the RSAE frequency sweeping process has also been observed in HL-2A
experiments Shi (2019).

In this work, we show that, as the RSAE frequency sweeps up during current
ramp, two important potential nonlinear mode coupling processes may happen. In
the first process, a RSAE
$\Omega_{R}\equiv\Omega_{R}\left(\omega_{R},\mathbf{k_{R}}\right)$ couples to
a TAE $\Omega_{T}\equiv\Omega_{T}\left(\omega_{T},\mathbf{k_{T}}\right)$ with
same poloidal and toroidal mode numbers and generates a geodesic acoustic mode
(GAM) $\Omega_{G}\equiv\Omega_{G}\left(\omega_{G},\mathbf{k_{G}}\right)$ with
toroidally symmetric and poloidally near symmetric mode structure Winsor et
al. (1968); Zonca and Chen (2008). Here, subscripts “R”, “T” and “G” represent
RSAE, TAE and GAM, respectively. In the second channel, a RSAE couples to a
counter-propagating TAE and generates an ion acoustic wave (IAW)
$\Omega_{S}\equiv\Omega_{S}\left(\omega_{S},\mathbf{k_{S}}\right)$, and this
process can occur even if RSAE has different poloidal/toroidal mode numbers
with TAE. Therefore, generally speaking, the condition for the latter process
to occur, including radial mode structure overlapping, is more easily
satisfied, as we show later. It is worthwhile noting that, in both channels,
the low frequency secondary modes, i.e., GAM and IAW, can both be Landau
damped due to resonance with thermal ions, and thus heat fuel deuterium and
tritium ions, proving a new alpha channelling mechanism that indirectly
transfer fusion alpha particle energy to fuel ions Fisch and Rax (1992); Hahm
(2015); Qiu et al. (2018).

We use the standard nonlinear perturbation theory to study the nonlinear
interactions between RSAE and TAE during current ramp. For typical low-$\beta$
discharges, magnetic compression can be neglected, and the electrostatic
potential $\delta\phi$ and the parallel component of vector potential $\delta
A_{\parallel}$ are used as perturbed field variables. Here, $\beta$ is the
ratio of thermal to magnetic pressure. Furthermore, we use
$\delta\psi\equiv\omega\delta A_{\parallel}/\left(ck_{\parallel}\right)$ as an
alternative field variable replacing $\delta A_{\parallel}$ for convenience,
and one can have $\delta\phi=\delta\psi$ in the ideal MHD limit.

For RSAE and TAE with $nq\gg 1$ in reactor relevant parameter regime Fasoli et
al. (2007); Chen and Zonca (2016); Wang et al. (2018); Ren et al. (2020), we
adopt the following ballooning mode representation in the
$\left(r,\theta,\varphi\right)$ field-aligned flux coordinates Connor et al.
(1978)

$\displaystyle\delta\phi$ $\displaystyle=$ $\displaystyle
Ae^{i\left(n\phi-\hat{m}\theta-\omega
t\right)}\underset{j}{\sum}e^{-ij\theta}\Phi\left(x-j\right)+c.c.{\color[rgb]{1,0,0},}{\color[rgb]{0,0,1}.}$

Here, $m=\hat{m}+j$ with $\hat{m}$ being the reference poloidal number,
$x\equiv nq-m$, $\Phi$ is the parallel mode structure with the typical radial
extension comparable to distance between neighboring mode rational surfaces,
and $A$ is the mode envelope amplitude. Furthermore, $\Omega_{G}$ and
$\Omega_{S}$ can be shown as

$\displaystyle\delta\phi_{G}$ $\displaystyle=$ $\displaystyle
A_{G}e^{i\left(\int k_{G}dr-\omega_{G}t\right)}+c.c.,$
$\displaystyle\delta\phi_{S}$ $\displaystyle=$ $\displaystyle
A_{S}e^{i\left(n_{S}\phi-m_{S}\theta-\omega_{S}t\right)}\Phi_{S}+c.c..$

The IAW is a secondary mode and its parallel mode structure $\Phi_{S}$ is
determined by $\Phi_{R}$ and $\Phi_{T}$ Qiu et al. (2016). Different from
previous works on spontaneously decay of a pump wave into two sidebands, e.g.,
as in the case of zonal flow excited by SAW Chen and Zonca (2012); Qiu et al.
(2017), the RSAE and TAE investigated here are simultaneously driven unstable
by EPs Wang et al. (2018), and thus, the radial mode structure of RSAE and TAE
may not overlap completely.

The nonlinear interaction of $\Omega_{R}$ and $\Omega_{T}$ can be described by
quasi-neutrality condition

$\displaystyle\dfrac{n_{0}e^{2}}{T_{i}}\left(1+\dfrac{T_{i}}{T_{e}}\right)\delta\phi_{k}$
$\displaystyle=$ $\displaystyle\underset{s}{\sum}\left\langle qJ_{k}\delta
H_{k}\right\rangle_{s},$ (1)

and nonlinear gyrokinetic vorticity equation derived from parallel Ampére’s
law Chen and Zonca (2016)

$\displaystyle\dfrac{c^{2}}{4\pi\omega_{k}^{2}}B\dfrac{\partial}{\partial
l}\dfrac{k_{\perp}^{2}}{B}\dfrac{\partial}{\partial
l}\delta\psi_{k}+\dfrac{e^{2}}{T_{i}}\left\langle\left(1-J_{k}^{2}\right)F_{0}\right\rangle\delta\phi_{k}$
$\displaystyle-\underset{s}{\sum}\left\langle\dfrac{q}{\omega_{k}}J_{k}\omega_{d}\delta
H_{k}\right\rangle$
$\displaystyle=-i\dfrac{c}{B_{0}\omega_{k}}\underset{\mathbf{k}=\mathbf{k^{\prime}}+\mathbf{k^{\prime\prime}}}{\sum}\mathbf{\hat{b}}\cdot\mathbf{k^{\prime\prime}}\times\mathbf{k^{\prime}}\left[\dfrac{c^{2}}{4\pi}k_{\perp}^{\prime\prime
2}\dfrac{\partial_{l}\delta\psi_{k^{\prime}}\partial_{l}\delta\psi_{k^{\prime\prime}}}{\omega_{k^{\prime}}\omega_{k^{\prime\prime}}}\right.$
$\displaystyle\left.+\left\langle
e\left(J_{k}J_{k^{\prime}}-J_{k^{\prime\prime}}\right)\delta
L_{k^{\prime}}\delta H_{k^{\prime\prime}}\right\rangle\right].$ (2)

Here, $J_{k}\equiv J_{0}\left(k_{\perp}\rho\right)$ with $J_{0}$ being the
Bessel function of zero index, $\rho=v_{\perp}/\Omega_{c}$ is the Larmor
radius with $\Omega_{c}$ being the cyclotron frequency, $F_{0}$ is the
equilibrium particle distribution function,
$\omega_{d}=\left(v_{\perp}^{2}+2v_{\parallel}^{2}\right)/\left(2\Omega_{c}R_{0}\right)\left(k_{r}\sin\theta+k_{\theta}\cos\theta\right)$
is the magnetic drift frequency, $l$ is the length along the equilibrium
magnetic field line, $\left\langle\cdots\right\rangle$ means velocity space
integration, $\sum_{s}$ is the summation of different particle species with
$s=i,e$ representing ion and electron, and $\delta
L_{k}\equiv\delta\phi_{k}-k_{\parallel}v_{\parallel}\delta\psi_{k}/\omega_{k}$.
The three terms on the left hand side of equation (2) are respectively the
field line bending, inertial and curvature coupling terms, dominating the
linear SAW physics. The two terms on the right hand side of equation (2)
correspond to Maxwell (MX) and Reynolds stresses (RS) that contribute to
nonlinear mode couplings as MX and RS doesn’t cancel each other Chen and Zonca
(2013), with their contribution dominating in the radially fast varying
inertial layer, and
$\sum_{\mathbf{k}=\mathbf{k^{\prime}}+\mathbf{k^{\prime\prime}}}$ indicates
the wavenumber and frequency matching condition required for nonlinear mode
coupling. $\delta H_{k}$ is the nonadiabatic particle response, which can be
derived from nonlinear gyrokinetic equation Frieman and Chen (1982):

$\displaystyle\left(-i\omega_{k}+v_{\parallel}\partial_{l}+i\omega_{d}\right)\delta
H_{k}=-i\omega_{k}\dfrac{q}{T}F_{0}J_{k}\delta L_{k}$
$\displaystyle-\dfrac{c}{B_{0}}\underset{\mathbf{k}=\mathbf{k^{\prime}}+\mathbf{k^{\prime\prime}}}{\sum}\mathbf{\hat{b}}\cdot\mathbf{k^{\prime\prime}}\times\mathbf{k^{\prime}}J_{k^{\prime}}\delta
L_{k^{\prime}}\delta H_{k^{\prime\prime}}.$ (3)

Here, with the diamagnetic drift related term neglected in the first term on
the right hand side of equation (3), we assume that thermal plasma -
dominating the nonlinear coupling process - contribution to RSAE/TAE linear
destabilization is negligible, and RSAE/TAE are excited by EPs. Inclusion of
thermal ion diamagnetic drift effect is straightforward, and will not change
the main physics picture here. For TAE/RSAE with
$\left|k_{\parallel}v_{e}\right|\gg\left|\omega_{k}\right|\gg\left|k_{\parallel}v_{i}\right|,\>\left|\omega_{d}\right|$,
the linear ion/electron responses can be derived to the leading order as
$\delta H^{L}_{k,i}=eF_{0}J_{k}\delta\phi_{k}/T_{i}$ and $\delta
H^{L}_{k,e}=-eF_{0}\delta\psi_{k}/T_{e}$. Furthermore, one can have, to the
leading order, ideal MHD constraint is satisfied, i.e.,
$\delta\phi_{T}=\delta\psi_{T}$, $\delta\phi_{R}=\delta\psi_{R}$, by
substituting these ion/electron responses of TAE and RSAE into quasi-
neutrality condition.

Figure 1: Cartoon of GAM generation by RSAE and TAE. The horizontal axis is
the normalized minor radius, the dashed purple curve is the $q$-profile, while
the red and blue curves are the $m/n=9/5$ and $m/n=10/5$ continua,
respectively. The reversed shear $q$-profile has a minimum value at $r\simeq
0.47$, where RSAE is generated above the local maximum of SAW continuum
induced by local $q_{min}$. TAE can be generated inside a radially nearby
continuum gap, with the frequency slightly higher than that of the RSAE. These
two modes can couple as they are radially overlapped, and generate a low
frequency GAM.

The nonlinear coupling and generation of GAM as the RSAE couples to the TAE
can be illustrated in Fig. 1, where a general $n=5$ SAW continuum in a typical
reversed shear configuration is given Bao et al. (2020). The $m=9$ and $m=10$
continuum are marked in red and blue, respectively, and they coupled at the
vicinity of $q=1.9$. A RSAE can be generated at the local maximum of the
continuum, while a TAE can be excited within the continuum gap, with their
radial localization determined by $q_{min}$ and $q=(2m+1)/(2n)$, respectively.
As $|nq_{min}-m|$ approaches $1/2$ during current ramp, the localization of
the RSAE gets closer to that of the TAE, and the frequency difference of RSAE
and TAE decreases. As a result, RSAE and TAE may couple, and generate a GAM.
Noting that GAM is characterized by the toroidal/poloidal mode numbers being
$n/m=0/0$ and a frequency much lower than those of RSAE/TAE, i.e.,
$|\omega_{G}|\ll|\omega_{R}|,|\omega_{T}|$, one can have that TAE and RSAE
have opposite poloidal and toroidal mode numbers as well as opposite
frequencies $\left(\omega_{R}\omega_{T}<0\right)$. Thus, the two modes
propagate in the same direction, i.e., their parallel phase velocities
$V_{P}\equiv\omega/k_{\parallel}$ have the same sign. For GAM with
predominantly electrostatic perturbation, its nonlinear generation can be
determined by the vorticity equation, while the quasi-neutrality condition is
used in determining the linear polarization of GAM and both AEs Qiu et al.
(2009); Cheng et al. (1985). The linear ion/electron responses to GAM can be
derived considering the
$\left|\omega_{tr,e}\right|\gg\left|\omega_{G}\right|\gg\left|\omega_{d,e}\right|$
and
$\left|\omega_{G}\right|\gg\left|\omega_{tr,i}\right|,\left|\omega_{d,i}\right|$
orderings, respectively, and one derives, to the leading order, $\delta
H_{G,i}^{L}=eF_{0}J_{G}\delta\phi_{G}/T_{i}$ and $\delta
H_{G,e}^{L}=-eF_{0}\overline{\delta\phi_{G}}/T_{e}$. Here, $\omega_{tr}\equiv
v_{\parallel}/(qR_{0})$ is the transit frequency, and
$\overline{(\cdots)}\equiv\int^{2\pi}_{0}(\cdots)d\theta/(2\pi)$ denotes flux
surface average. The nonlinear equation describing GAM generation can then be
derived from nonlinear vorticity equation as

$\displaystyle\dfrac{e^{2}}{T_{i}}\left\langle\left(1-J_{G}^{2}\right)F_{0}\right\rangle\delta\phi_{G}-\underset{s}{\sum}\left\langle\dfrac{q}{\omega}J_{G}\omega_{d}\delta
H^{L}_{G}\right\rangle$
$\displaystyle\simeq-i\dfrac{c}{B_{0}\omega}\mathbf{\hat{b}}\cdot\mathbf{k_{R}}\times\mathbf{k_{T}}$
$\displaystyle\times\left[\dfrac{c^{2}}{4\pi}\left(k_{T\perp}^{2}-k_{R\perp}^{2}\right)\dfrac{k_{T,\parallel}k_{R,\parallel}}{\omega_{T}\omega_{R}}\delta\psi_{T}\delta\psi_{R}\right.$
$\displaystyle\left.+\left\langle
e\left(J_{T}-J_{R}\right)\left(\delta\phi_{T}\delta H_{R}+\delta\phi_{R}\delta
H_{T}\right)\right\rangle\right].$ (4)

In deriving equation (4), the linear field line bending term associated with
the GAM electromagnetic effects is neglected because GAM is predominantly an
electrostatic mode. The last term of equation (4) is the RS evaluated using
the linear ion responses to RSAE/TAE and the $k^{2}_{\perp}\rho^{2}_{i}\ll 1$
ordering. Noting that $k_{G}=k_{T,r}+k_{R,r}$, $k_{r}\gg k_{\theta}$ for
inertial layer responses of TAE and RSAE that dominates the nonlinear mode
coupling Chen and Zonca (2013), one obtains

$\displaystyle\delta\phi_{G}$ $\displaystyle=$
$\displaystyle\dfrac{c\omega}{B_{0}\left(\omega^{2}-\omega_{G}^{2}\right)}k_{T,\theta}\left(1-\dfrac{k_{T,\parallel}k_{R,\parallel}}{\omega_{T}\omega_{R}}v_{A}^{2}\right)$
(5) $\displaystyle\times$
$\displaystyle\left(\delta\phi_{R}\partial_{r}\delta\phi_{T}-\delta\phi_{T}\partial_{r}\delta\phi_{R}\right).$

Here, $\omega_{G}$ is the eigenfrequency of GAM and is given as
$\omega^{2}_{G}=\left(7/4+T_{e}/T_{i}\right)v_{ti}^{2}/R_{0}^{2}$ Qiu et al.
(2009), $\omega=\omega_{T}+\omega_{R}$ from balancing the temporal evolution
of both sides of equation (4), which is not necessarily exactly the same as
$\omega_{G}$. Equation (5) describes the nonlinear drive of GAM by the beating
of co-existing TAE and RSAE simultaneously driven unstable by e.g., EPs. The
effective coupling and GAM generation, requires that
$1-k_{T,\parallel}k_{R,\parallel}v^{2}_{A}/(\omega_{T}\omega_{R})\neq 0$,
i.e., breaking of pure Alfvénic state Chen and Zonca (2013). This condition is
satisfied by the deviation of both RSAE and TAE from ideal SAW dispersion
relation, due to the effects of reversed shear profile and toroidicity Zonca
et al. (2002); Zonca and Chen (1996), respectively, and one has
$1-k_{T,\parallel}k_{R,\parallel}v^{2}_{A}/(\omega_{T}\omega_{R})\simeq
O(\epsilon)$ with $\epsilon\equiv r/R_{0}$ being the inverse aspect ratio.
Noting that, this forced driven process is thresholdless, but only when
$\omega_{T}+\omega_{R}$ gets sufficiently close to $\omega_{G}$ during RSAE
frequency sweeping as a result of current ramp, GAM will be strongly excited,
and the process can be observed experimentally and play an important role in
the RSAE/TAE saturation.

In Ref. Van Zeeland et al. (2007), a central localized RSAE coupling to TAE
and generate a global TAE is investigated, where a continuum similar to Fig. 1
is used for the illustration of the mechanism. However, the process discussed
in Ref. Van Zeeland et al. (2007) is a linear process, where the RSAE has the
same toroidal mode number with the TAE while neighboring poloidal mode number
is coupled through toroidicity, and is thus completely different from the
present work.

Figure 2: Cartoon of IAW generation by RSAE and TAE. The dashed purple curve
is the $q$-profile, while the red and blue curves are SAW continua of $n=5$
and $n=7$. An $n=5$ RSAE localizes at the local maximum of the SAW continuum
caused by the shear reversal, and an $n=7$ TAE localizes within the SAW
continuum induced by toroidicity.

Considering the other case of RSAE coupling to TAE and generating an IAW with
finite $k_{S,\parallel}$ and $|\omega_{S}|\ll|\omega_{T}|,|\omega_{R}|$, one
can have TAE and RSAE have opposite frequencies while their parallel
wavenumbers have the same sign, and thus are counter-propagating with opposite
$V_{p}$. This process is illustrated in Fig. 2 using the nonlinear coupling of
an $n=5$ RSAE and an $n=7$ TAE as example, where the red and blue curve
represent the $n=5$ and $n=7$ SAW continua, and the corresponding RSAE and TAE
are radially overlapped, leading to potential nonlinear coupling. Note that,
the RSAE and TAE can have different toroidal mode numbers, unlike the case of
GAM generation limited by the $n_{R}+n_{T}=0$ constraint. Thus, it is expected
the condition for IAW generation can be more easily satisfied (as shown
clearly in Fig. 2), in future burning plasmas where multiple-$n$ SAW
instabilities can be driven unstable simultaneously.

Separating the linear to nonlinear particle responses to IAW as $\delta
H_{S}=\delta H_{S}^{L}+\delta H_{S}^{NL}$, the leading order linear
ion/electron responses to IAW can be derived from gyrokinetic equation as
$\delta
H_{S,i}^{L}=\left(eF_{0}J_{S}\omega_{S}\delta\phi_{S}\right)/\left[T_{i}\left(\omega_{S}-k_{S,\parallel}v_{\parallel}\right)\right]$
and $\delta H_{S,e}^{L}=0$, considering the
$\left|\omega_{tr,e}\right|\gg\left|\omega_{S}\right|\simeq\left|\omega_{tr,i}\right|\gg\left|\omega_{d,i}\right|,\left|\omega_{d,e}\right|$
ordering. Furthermore, the nonlinear particle responses to IAW can be derived
as

$\displaystyle\delta H_{S,e}^{NL}$ $\displaystyle=$
$\displaystyle-i\dfrac{\hat{\Lambda}}{\omega_{T}}\dfrac{e}{T_{e}}F_{0}\delta\psi_{T}\delta\psi_{R},$
(6) $\displaystyle\delta H_{S,i}^{NL}$ $\displaystyle=$
$\displaystyle-i\dfrac{\hat{\Lambda}}{\omega_{T}}\dfrac{e}{T_{i}}F_{0}J_{T}J_{R}\delta\psi_{T}\delta\psi_{R},$
(7)

with
$\hat{\Lambda}=\left(c/B_{0}\right)\mathbf{\hat{b}}\cdot\mathbf{k_{R}}\times\mathbf{k_{T}}$.
Here, $\omega_{S}\simeq-\omega_{T}$ is used. Assuming ideal MHD condition
($\delta\psi=\delta\phi$) for RSAE and TAE in the radially fast varying
inertial layer where nonlinear coupling dominates, the excitation of
electrostatic IAW by TAE and RSAE can be derived from quasi-neutrality
condition, by substituting the particle responses into quasi-neutrality
condition, and one obtains

$\displaystyle\mathscr{E}_{S}\delta\phi_{S}$ $\displaystyle=$ $\displaystyle
i\dfrac{\hat{\Lambda}}{\omega_{T}}\beta_{S}\delta\phi_{T}\delta\phi_{R}.$ (8)

Here, $\mathscr{E}_{S}\equiv 1+\tau+\tau\left\langle
F_{0}J_{S}^{2}/n_{0}\right\rangle\zeta_{S}Z\left(\zeta_{S}\right)$ is the
linear dispersion function of $\Omega_{S}$, with $\tau\equiv T_{e}/T_{i}$,
$\zeta_{S}\equiv\omega_{S}/\left(k_{S,\parallel}v_{i}\right)$ and
$Z\left(\zeta_{S}\right)$ being the plasma dispersion function defined as

$\displaystyle Z\left(\zeta_{S}\right)$
$\displaystyle\equiv\frac{1}{\sqrt{\pi}}\int_{-\infty}^{\infty}\frac{e^{-y^{2}}}{y-\zeta_{S}}dy.$

In addition, $\beta_{S}\equiv 1+\tau\left\langle
F_{0}J_{S}J_{T}J_{R}/n_{0}\right\rangle\left(1+\zeta
Z\left(\zeta\right)\right)$ is related to the cross-section of the coupling,
and in the case of counter-propagating RSAE/TAE of interest here, the RS and
MX have the same sign, and $\beta_{S}\simeq 2$ in the long wavelength limit.
Similarly with the case discussed before, equation (8) describes the nonlinear
generation of IAW by the coupling of TAE and RSAE, showing that this coupling
process is thresholdless, but only significantly affect the nonlinear dynamics
as $|\omega_{R}+\omega_{T}|\lesssim O(v_{i}/(qR_{0}))$ (i.e.,
$\zeta_{S}\lesssim O(1)$). In this parameter regime, RSAE and TAE strongly
couples as $|\mathscr{E}_{S}|\ll 1$, as shown by equation (8). On the other
hand, as $\zeta_{S}\lesssim O(1)$, the nonlinearly generated IAW is a heavily
Landau damped quasi-mode that effectively gives energy to thermal ions,
leading to effective dissipation of RSAE/TAE wave energy.

In conclusion, two novel nonlinear coupling channels for RSAE saturation are
proposed and investigated. In the first channel, the RSAE couples to a co-
propagating TAE with the same toroidal and poloidal mode numbers, and
generates a GAM. In the second channel, the RSAE couples to a counter-
propagating TAE and generates an IAW. In the latter case, the RSAE and TAE
toroidal/poloidal mode numbers are not necessarily the same. The frequency
matching condition is more easily satisfied as the RSAE frequency is slightly
lower than that of TAE, and is more easily satisfied during current ramp. Our
results show that, the IAW generation can be more efficient in that 1. the
RSAE and TAE don’t have to have the same toroidal/poloidal mode numbers, and
thus, the condition for mode structure radial overlapping can be more easily
satisfied, and 2. the RSAE and TAE are counter-propagating, so the RS and MX
will not cancel each other, and thus, the nonlinear coupling cross-section is
much larger. In burning plasmas, the RSAE and TAE are linearly excited by EPs
simultaneously, and as a result, both processes are thresholdless. In both
processes, the nonlinearly generated secondary mode, i.e., GAM and IAW, can be
Landau damped due to resonance with thermal ions, and heat thermal ions, and
thus provides an effectively channel for nonlinearly transfer fusion alpha
particles to fuel ions for the sustained burning of fusion plasmas. The
saturation level of the low frequency sideband in the two processes, can be
estimated from equation(5) or (8), by expanding the GAM or IAW linear
dispersion function along the characteristics Qiu et al. (2018) and balancing
the nonlinear drive by RSAE&TAE nonlinear coupling, and the GAM/IAW Landau
damping. The corresponding fuel ion anomalous heating rate can also be
obtained, and will be investigated in a future publication.

This work is supported by the National Key R&D Program of China under Grant
No. 2017YFE0301900, the National Science Foundation of China under grant No.
11875233, and China Postdoctoral Science Foundation under grant No.
2020M670756. The authors are grateful to Dr. Jian Bao (Institute of physics,
Chinese Academy of Science) for his help with the SAW continua.

## References

  * Tomabechi et al. (1991) K. Tomabechi, J. Gilleland, Y. Sokolov, R. Toschi, and I. Team, Nuclear Fusion 31, 1135 (1991).
  * Fasoli et al. (2007) A. Fasoli, C. Gormenzano, H. Berk, B. Breizman, S. Briguglio, D. Darrow, N. Gorelenkov, W. Heidbrink, A. Jaun, S. Konovalov, et al., Nuclear Fusion 47, S264 (2007).
  * Chen and Zonca (2016) L. Chen and F. Zonca, Review of Modern Physics 88, 015008 (2016).
  * Chen and Wang (2020) W. Chen and Z. Wang, Chinese Physical Letters 37, 125001 (2020).
  * Ding et al. (2015) R. Ding, R. Pitts, D. Borodin, S. Carpentier, F. Ding, X. Gong, H. Guo, A. Kirschner, M. Kocan, J. Li, et al., Nuclear Fusion 55, 023013 (2015).
  * Chen (1994) L. Chen, Physics of Plasmas 1, 1519 (1994).
  * Cheng et al. (1985) C. Cheng, L. Chen, and M. Chance, Ann. Phys. 161, 21 (1985).
  * Chen (1988) L. Chen, in _Theory of Fusion Plasmas_ , edited by J. Vaclavik, F. Troyon, and E. Sindoni (Association EUROATOM, Bologna, 1988), p. 327.
  * Fu and Van Dam (1989) G. Y. Fu and J. W. Van Dam, Physics of Fluids B 1, 1949 (1989).
  * Berk et al. (2001) H. L. Berk, D. N. Borba, B. N. Breizman, S. D. Pinches, and S. E. Sharapov, Phys. Rev. Lett. 87, 185002 (2001).
  * Zonca et al. (2002) F. Zonca, S. Briguglio, L. Chen, S. Dettrick, G. Fogaccia, D. Testa, and G. Vlad, Physics of Plasmas 9, 4939 (2002).
  * Kramer et al. (2004) G. J. Kramer, N. N. Gorelenkov, R. Nazikian, and C. Z. Cheng, Plasma Physics and Controlled Fusion 46, L23 (2004).
  * Heidbrink et al. (1993) W. Heidbrink, E. Strait, M. Chu, and A. Turnbull, Phys. Rev. Lett. 71, 855 (1993).
  * Buratti et al. (2005) P. Buratti, P. Smeulders, F. Zonca, S. Annibaldi, M. De Benedetti, H. Kroegler, G. Regnoli, O. Tudisco, et al., Nuclear Fusion 45, 1446 (2005).
  * Hahm and Chen (1995) T. S. Hahm and L. Chen, Phys. Rev. Lett. 74, 266 (1995).
  * Zonca et al. (1995) F. Zonca, F. Romanelli, G. Vlad, and C. Kar, Phys. Rev. Lett. 74, 698 (1995).
  * Todo et al. (2010) Y. Todo, H. Berk, and B. Breizman, Nuclear Fusion 50, 084016 (2010).
  * Biancalani et al. (2010) A. Biancalani, L. Chen, F. Pegoraro, and F. Zonca, Physical Review Letters 105, 095002 (2010).
  * Chen and Zonca (2012) L. Chen and F. Zonca, Phys. Rev. Lett. 109, 145002 (2012).
  * Zhang and Lin (2013) H. Zhang and Z. Lin, Plasma Science and Technology 15, 969 (2013).
  * Biancalani et al. (2018) A. Biancalani, A. Bottino, A. D. Siena, T. Gorler, I. Novikau, O. Gurcan, P. Morel, F. Jenko, E. Lanti, N. Ohana, et al. (Ahmedabad, India, 2018).
  * Gong et al. (2019) X. Gong, A. Garofalo, J. Huang, J. Qian, C. Holcomb, A. Ekedah, R. Maingi, E. Li, L. Zeng, B. Zhang, et al., Nuclear Fusion 59, 086030 (2019).
  * Huang et al. (2020) J. Huang, A. Garofalo, J. Qian, X. Gong, S. Ding, J. Varela, J. Chen, W. Guo, K. Li, M. Wu, et al., Nuclear Fusion 60, 126007 (2020).
  * Wang et al. (2018) T. Wang, Z. Qiu, F. Zonca, S. Briguglio, G. Fogaccia, G. Vlad, and X. Wang, Physics of Plasmas 25, 062509 (2018).
  * Ren et al. (2020) Z. Ren, Y. Chen, G. Fu, and Z. Wang, Nuclear Fusion 60, 016009 (2020).
  * Van Zeeland et al. (2007) M. A. Van Zeeland, M. E. Austin, N. N. Gorelenkov, W. W. Heidbrink, G. J. Kramer, M. A. Makowski, G. R. McKee, R. Nazikian, E. Ruskov, and A. D. Turnbull, Physics of Plasmas 14, 056102 (2007).
  * Wang et al. (2020) T. Wang, Z. Qiu, F. Zonca, S. Briguglio, and G. Vlad, Nuclear Fusion 60, 126032 (2020).
  * Shi (2019) P. Shi (2019), private communication.
  * Winsor et al. (1968) N. Winsor, J. L. Johnson, and J. M. Dawson, Physics of Fluids 11, 2448 (1968).
  * Zonca and Chen (2008) F. Zonca and L. Chen, Europhys. Lett. 83, 35001 (2008).
  * Fisch and Rax (1992) N. J. Fisch and J.-M. Rax, Phys. Rev. Lett. 69, 612 (1992).
  * Hahm (2015) T. S. Hahm, Plasma Science and Technology 17, 534 (2015).
  * Qiu et al. (2018) Z. Qiu, L. Chen, F. Zonca, and W. Chen, Phys. Rev. Lett. 120, 135001 (2018).
  * Connor et al. (1978) J. Connor, R. Hastie, and J. Taylor, Phys. Rev. Lett. 40, 396 (1978).
  * Qiu et al. (2016) Z. Qiu, L. Chen, and F. Zonca, Nuclear Fusion 56, 106013 (2016).
  * Qiu et al. (2017) Z. Qiu, L. Chen, and F. Zonca, Nuclear Fusion 57, 056017 (2017).
  * Chen and Zonca (2013) L. Chen and F. Zonca, Physics of Plasmas 20, 055402 (2013).
  * Frieman and Chen (1982) E. A. Frieman and L. Chen, Physics of Fluids 25, 502 (1982).
  * Bao et al. (2020) J. Bao, W. Zhang, D. Li, and Z. Lin, _Effects of plasma diamagnetic drift on alfvén continua and discrete eigenmodes in tokamaks_ (2020).
  * Qiu et al. (2009) Z. Qiu, L. Chen, and F. Zonca, Plasma Physics and Controlled Fusion 51, 012001 (2009).
  * Zonca and Chen (1996) F. Zonca and L. Chen, Physics of Plasmas 3, 323 (1996).


# A Response to Glaze Purification via IMPRESS

Shawn Shan†, Stanley Wu†, Haitao Zheng, Ben Y. Zhao
† denotes authors with equal contribution
Department of Computer Science, University of Chicago
{shawnshan, stanleywu, htzheng<EMAIL_ADDRESS>
###### Abstract.

Recent work proposed a new mechanism to remove protective perturbation added
by Glaze in order to again enable mimicry of art styles from images protected
by Glaze. Despite promising results shown in the original paper, our own tests
with the authors’ code demonstrated several limitations of the proposed
purification approach. The main limitations are 1) purification has a limited
effect when tested on artists that are not well-known historical artists
already embedded in original training data, 2) problems in evaluation metrics,
and 3) collateral damage on mimicry result for clean images. We believe these
limitations should be carefully considered in order to understand real world
usability of the purification attack.

††copyright: none

##  1\. Introduction

Cao et al.  (Cao et al., 2023) recently proposed a purification-based method
(IMPRESS) to remove Glaze protection (Shan et al., 2023) that are designed to
protect art styles from AI mimicry. We first summarize how Glaze and IMPRESS
work in §2. Then we show IMPRESS achieve limited effectiveness against Glaze
in a variety of real-world mimicry scenarios. In short:

  * •

IMPRESS has weaker performance on artists who are not historical artists
already with a large presence in the pretrained models;

  * •

IMPRESS works poorly on smooth surface art styles (e.g., realism);

  * •

IMPRESS damages image quality even for clean images.

##  2\. Background

Here, we first summarize Glaze and then IMPRESS purification method.

Glaze protection against style mimicry.  Glaze seeks to protect artist’s
artwork from AI mimicry by adding small perturbations on these artwork to
confuse diffusion models. Given an artwork $x$ and target style $T$ that is
different from the artist’s, Glaze first uses a pretrained style transfer
model $\Omega$ to compute a style transferred version of the artwork. We
denote such image as $\Omega(x,T)$. Then, Glaze computes a cloak $\delta_{x}$
that optimize the latent representation of Glazed artwork ($x+\delta_{x}$) to
be similar to the style transferred artwork ($\Omega(x,T)$). The Glaze
optimization effectively moves the original image to a new position in the
high dimensional latent space, causing model to learn an incorrect art style.
Glaze calculates the latent space using the feature extractor ($\mathcal{E}$)
from a diffusion model. Formally, we write the Glaze optimization as solving
the following:

(1)
$\displaystyle\min_{\delta_{x}}||\mathcal{E}(x+\delta_{x})-\mathcal{E}(\Omega(x,T))||_{2}$
$\displaystyle\text{subject to }\text{LPIPS}(x+\delta_{x},x)<p_{G}$

We use LPIPS, a popular human-perceived visual distortion metric (Zhang et
al., 2018), to bound the perturbation within a budget $p_{G}$.

IMPRESS Purification Method.  IMPRESS adds additional perturbation on top of a
Glazed artwork hoping to “purify” the Glaze effect – recovering the precise
latent representation of original artwork. First, the authors empirically find
that when passing Glazed images through an image autoencoder, the output image
looks more different from the input image, compared to the output when
inputting a clean image to the same autoencoder. Then authors assume removing
this particular discrepancy would guide them to find the original (non-Glazed)
image.

IMPRESS purification optimizes perturbations on Glaze images such that
purified images behave similarly to clean images when passing through an
autoencoder. The authors assume the optimization process will guide the image
to move back to the original latent space of the non-Glazed image. Formally,
IMPRESS purification optimize a perturbation $\delta_{pur}$ on a Glaze image
$x_{glazed}$:

(2)
$\displaystyle\min_{\delta_{pur}}||(x_{glazed}+\delta_{pur})-\text{VAE}(x_{glazed}+\delta_{pur})||_{2}^{2}$
$\displaystyle\text{subject to
}\text{LPIPS}(x_{glazed}+\delta_{pur},x_{glazed})<p_{I}$

VAE is an image autoencoder, which consists of an encoder $\mathcal{E}$
followed by a decoder $D$. IMPRESS uses the same autoencoder as the stable
diffusion model. The perturbation $\delta_{pur}$ is bounded by a LPIPS
perturbation budget $p_{I}$.

##  3\. Evaluation of IMPRESS

Here, we first identify a few flaws in IMPRESS’s experiment setup and then
evaluate purification performance in several realistic mimicry scenarios.

Setup.  We follow the purification setup from the original IMPRESS paper and
the source code provided by the authors. We use default Glaze setup
(perturbation budget = $0.07$) as described in the paper. We follow style
mimicry setup in prior work (Shan et al., 2023). Given a set of artworks, we
mimic its style by fine-tuning stable diffusion model (version 1.5) using
DreamBooth for 600 to 1000 steps depending on the finetuning sample size for
the artist.

Figure 1. Mimicry results on non-historical artists (Karla Ortiz and Kim van
Deun). Left shows the images generated from a model trained on original
images; right shows the images generated from a model trained on images that
are first Glazed and then purified by IMPRESS.  Figure 2. Mimicry results on
smooth surface art styles (Ivan Shishkin and Ivan Aivazovsky). Left shows the
images generated from a model trained on original images; right shows the
images generated from a model trained on images that are first Glazed and then
purified by IMPRESS.  Figure 3. Original artwork and corresponding purified
artwork.

###  3.1. Generalizability to Realistic Scenarios

In the original IMPRESS paper, the authors focus the evaluation on protecting
the art styles of famous historical artists (Monet, Van Gogh) – whose style
are already learned by pretrained diffusion models prior to style mimicry. In
the real-world, it is current artists who are most concerned about AI
mimicking their art style. Glaze is designed to protect those artists, and
they are not as heavily pretrained into the base model as Monet or Van Gogh.

We evaluate IMPRESS on non-historical artists and show that purification has
limited effectiveness at removing Glaze protection. Even for historical
artists, we find the purification works poorly for certain art styles. Lastly,
we find purification also degrades clean image quality where it removes
texture from art pieces.

Performance on non-historical artists.  We use artwork from Karla Ortiz and
Kim Van Deun to simulate the mimicry attack on current artists. Karla is a
fine-art artist with a similar style as some historical artists tested in the
IMPRESS paper. Figure 1 shows the mimicry results when model trained on
artist’s original art pieces (left) or when model trained on art pieces that
are first Glazed and then purified by IMPRESS (right). We can observe
significant amount of artifacts on mimicry results when the model is trained
on purified Glazed images.

Poor performance on certain art styles.  IMPRESS works by adding artifacts on
Glazed images to recover the latent representation of original artwork. We
found purification has more challenges recovering smooth surface art styles
(realism art, symbolism, romanticism, etc) even for historical artists already
trained into the base model. We choose two historical artists: Ivan Aivazovsky
(romanticism style) and Ivan Shishkin (realism style).

Figure 2 shows the mimicry results. We see IMPRESS introduces signficiant
amount of artifacts to the mimicked images. The weaker performance is likely
because purifying the original smoother surface art requires the optimization
to be very percise – find the exact smooth surface.

Degrading image quality.  We found the purification process degrades clean
image quality. Figure 3 shows original artwork and its corresponding purified
artwork. The purified artwork is much more blurry as purification removes
textures from the images.

###  3.2. CLIP-based metrics are Inaccurate

“CLIP genre accuracy” quantifies if mimicked art is classified into the same
art genre as the original art pieces according to a CLIP model. It has been
used in prior work to evaluate mimicry success. However, in our own tests
dating back to late September 2023, we found CLIP genre accuracy is a poor
indicator of end-to-end mimicry success. CLIP accuracy is especially poor when
evaluating attacks against Glaze. The reason is that attacks (such as IMPRESS
purification), as they seek to remove Glaze effect from art, often have
signficant impact on the base image quaility of the artwork. The degradation
in image quality is not captured by CLIP accuracy, e.g., a very low quality
cubism painting is still classified as “cubism” with high probability. But the
result are not successful mimicry due to the low quaility of the mimicked
images. Because of these poor properties as an accuracy indicator, we stopped
using CLIP distance as a success metric starting with the Glaze v1.1 update
(October 2023).

##  4\. Discussion

We conclude with a discussion of Glaze, then describe some potential
implications of attempts to remove Glaze.

Challenges of evaluating mimicry protection.  Evaluating mimicry performance
is challenging with automated metrics. As we shown in §3, CLIP-based metric do
not work well. Another popular image quaility metric, Fréchet inception
distance (FID), is also shown to be erronous at evaluating generative models
(Podell et al., 2023). The standard evaluation approach today is to rely on
human inspectors (Podell et al., 2023), which is not only expensive, but also
produces results with high variance. Particularly for art mimicry, the average
human user (found on MTurk or Profilic) often does not have sufficient
expertise to judge whether a mimicry is successful or not, while sourcing to
professional artists can be prohibitively expensive. We hope future work can
address this challenge by designing automated metrics that can faithfully
approximate ratings from human artists.

Implications on Copyright Law. Technical research results and their
applications rarely have legal implications in practice. In the case of Glaze
and image copyrights, however, there are interesting open questions with
respect to copyright law. More specifically, Section 1201 of US Copyright law
(enacted in 1998 as part of the DMCA) “prohibits circumvention of
technological protection measures employed by or on behalf of copyright owners
to control access to their works.”111“Legal Liability for Removal or
Alteration of Copyright Management Information Under the DMCA,” John
DiGiacomo, September 2020, https://revisionlegal.com/copyright/legal-
liability-for-removal-or-alteration-of-copyright-management-information-under-
the-dmca/.

For artists who apply Glaze to prevent unauthorized use of their art images
for training, applying Glaze to one’s art likely qualifies as an explicit
statement of “opt-out” and a measure to control access to copyrighted work. In
this uncertain legal and regulatory landscape, it is an interesting legal
question to explore whether Section 1201 applies to the use of Glaze and
attempted counter-measures to Glaze.

## Acknowledgements

We are grateful to Bochuan Cao, Changjiang Li, Ting Wang, Jinyuan Jia, Bo Li,
Jinghui Chen for providing us code and discussing their attack with us prior
to paper publication.

## References

  * (1)
  * Cao et al. (2023) Bochuan Cao, Changjiang Li, Ting Wang, Jinyuan Jia, Bo Li, and Jinghui Chen. 2023\.  IMPRESS: Evaluating the Resilience of Imperceptible Perturbations Against Unauthorized Data Usage in Diffusion-Based Generative AI.  _Proc. of NeurIPS_.
  * Podell et al. (2023) Dustin Podell, Zion English, Kyle Lacey, Andreas Blattmann, Tim Dockhorn, Jonas Müller, Joe Penna, and Robin Rombach. 2023\.  Sdxl: Improving latent diffusion models for high-resolution image synthesis.  _arXiv preprint arXiv:2307.01952_ (2023).
  * Shan et al. (2023) Shawn Shan, Jenna Cryan, Emily Wenger, Haitao Zheng, Rana Hanocka, and Ben Y Zhao. 2023\.  Glaze: Protecting artists from style mimicry by text-to-image models.  _Proc. of USENIX Security_.
  * Zhang et al. (2018) Richard Zhang, Phillip Isola, Alexei A Efros, Eli Shechtman, and Oliver Wang. 2018.  The unreasonable effectiveness of deep features as a perceptual metric. In _Proc. of CVPR_. 586–595.


# Data-Driven Predictive Control for Continuous-Time Industrial Processes with
Completely Unknown Dynamics

Yuanqiang Zhou, Dewei Li, Yugeng Xi Y. Zhou, D. Li, Y. G. Xi are with the
Department of Automation, Shanghai Jiao Tong University and the Key Laboratory
of System Control and Information Processing, Ministry of Education of China,
Shanghai, 200240, China, email<EMAIL_ADDRESS>work was supported by National Natural Science Foundation of China (NSFC)
under Grant No. 61333009, 61521063, 61473317, 61573239. The work of Y. Zhou is
supported by the China Scholarship Council.

###### Abstract

This paper investigates the data-driven predictive control problems for a
class of continuous-time industrial processes with completely unknown
dynamics. The proposed approach employs the data-driven technique to get the
system matrices online, using input-output measurements. Then, a model-free
predictive control approach is designed to implement the receding-horizon
optimization and realize the reference tracking. Feasibility of the proposed
algorithm and stability of the closed-loop control systems are analyzed,
respectively. Finally, a simulation example is provided to demonstrate the
effectiveness of the proposed approach.

###### Index Terms:

Data-driven control, industrial processes, model predictive control (MPC),
reference tracking.

##  I Introduction

As a practically effective approach, model predictive control (MPC), or
receding horizon control, has attracted notable attention in the field of
industrial process control [1]. To deal with optimal control problems, MPC can
allow for industrial processes uncertainties and constraints much more
straightforwardly than other methods [2, 3]. The core of all model-based
predictive algorithms is to use “open-loop optimal control” instead of
“closed-loop optimal control” within a moving horizon [4]. It brings a lot of
robustness and reliability to allow the controller have the ability to
recognize the control process. But, the acquisition of knowledge of a priori
model seriously affects those performances of MPC. The practically measured
data comes from the complicated processes and its utilization into MPC will
greatly facilitate the design procedures, avoiding the need for initially
accurate dynamic models [5].

Nowadays, many research works focus on the data-driven predictive control,
that is applicable on-line both to regulation and tracking control problems
[6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]. More earlier works can refer to
[18] and [13]. A Markov data-based LQG control algorithm is suggested in [18]
and data-based optimal control based on the system’s Markov parameters is
provided in [13]. Both of them utilize the prior measurements to design the
predictive control, along with its implementation on-line. Many results (see,
[9, 10, 11]) still need some knowledge of the system information to fit the
structure of the process model. Other works, like [19], focus on the polytopic
uncertain systems with unmeasurable system states and a convex hull needs to
be known. In [16, 20], data-driven subspace approach is introduced to design
the predictive controller and in [21, 8], reinforcement learning approach is
used to reduce the model-based dependence on predictive controller design
procedures. We remark that all the aforementioned predictive control methods
are designed based on the accurate model, as well as adequate uncertainty
description of the linear or non-linear plant of the processes.

In this paper, the limitations of them are circumvented. We use the adequately
measured data from the complicated industrial processes following the methods
presented in [22, 23] and [24, 25]. To be precise, the data-driven learning
technique of [22] will be employed to iteratively approximate the dynamical
parameters, without requiring the _a prior_ knowledge of the system matrices.
Then, the linear plant’s version of [26] will be applied to predict the future
trajectories, by following the continuous-time predictive control approach of
[7] but removing the assumption on partial knowledge of the system dynamics.
Under this framework, the data-driven predictive control input can be
generated on-line and can be used for the control of time-varying or nonlinear
plants, since the algorithm is able to adapt to the actual dynamics by
obtaining a linear model of the system at each sample. The contributions of
the paper are three-fold. (i) A data-driven approach is proposed to adaptively
approximate the system matrices, without requiring the _a prior_ knowledge of
the system. (ii) A continuous-time data-driven MPC approach is developed for
the continuous-time linear system, using repeatedly the state and input
formation on some fixed time intervals. (iii) By implementing the proposed
data-driven predictive control algorithm, both recursive feasibility of the
optimization problem and closed-loop stability of the whole system are
guaranteed.

The rest of this paper is organized as follows. In Section II, the problems
are briefly formulated. In Section III, the data-driven predictive control
approach with completely unknown dynamics algorithm is presented, and the
feasibility and stability analysis are conducted. In Section IV, we apply the
proposed approach to the optimal control problem of two continuous stirred
tank reactor (CSTR). Conclusions are given in Section V.

Notation: Through this note, $\mathbb{R}$ denotes the set of real numbers,
$|\cdot|$ represents the Euclidean norm for a vector and the induced norm for
a matrix. For a $a\in\mathbb{R}$, $\mathbb{R}_{\geq a}$ denotes the interval
$[a,\infty)$ and $\mathbb{R}_{>a}$ denotes the interval $(a,\infty)$. For any
vector $x$, $x^{\textrm{T}}$ denotes its transpose and ${\|x\|^{2}_{P}}$
represents ${{x^{\textrm{T}}}Px}$ for a real symmetrix and positive defined
matrix $P$. For a real symmetric matrix $A\in\mathbb{R}^{n\times n}$, $A\succ
0$ or $A\succeq 0$ means that $A$ is positive definite or semi-positive
definite, $\lambda_{M}(A)$ and $\lambda_{m}(A)$ denote the maximum and minimum
eigenvalue of $A$, respectively. $\otimes$ indicates the Kronecker product
operator, $\text{vec}(\cdot)$ denotes the vectorization operator and
$\text{vec}^{-1}(\cdot)$ denotes the converse vectorization operator, i.e.,
$\text{vec}(A)=[a_{1}^{\textrm{T}},\ldots,a_{m}^{\textrm{T}}]$ and
$\text{vec}^{-1}[a_{1}^{\textrm{T}},\ldots,a_{m}^{\textrm{T}}]=A$, where
$a_{i}\in\mathbb{R}^{n}$ are the columns of a matrix $A\in\mathbb{R}^{n\times
m}$. A continuous function $\alpha:\mathbb{R}_{\geq
0}\rightarrow\mathbb{R}_{\geq 0}$ is said to be a $\mathcal{K}$ function if it
is strictly increasing, and $\alpha(s)>0$ for $s>0$ with $\alpha(0)=0$. A
continuous function $\alpha(\cdot)$ is said to be a $\mathcal{K}_{\infty}$
function if it is a $\mathcal{K}$ function, and $\alpha(s)\rightarrow\infty$
for $s\rightarrow 0$. For any piecewise continuous function
$u(\cdot):\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}^{m}$, $u^{[\textrm{r}]}$
denotes the $r$th order derivatives of $u(\cdot)$.

##  II Problem formulation and preliminaries

Consider a continuous-time industrial process control system with the form

$\displaystyle{\dot{x}}(t)=$ $\displaystyle A{x}(t)+Bu(t),\quad
x(t_{0})=x_{0},$ (1) $\displaystyle y(t)=$ $\displaystyle Cx(t),$ (2)

where $t\geq{t_{0}}$, ${x}\in\mathbb{R}^{n}$ is the measurable state;
$u\in\mathcal{U}\subset\mathbb{R}^{m}$ is the control input;
$y\in\mathbb{R}^{q}$ is the controlled output fully available for feedback
control design. Assume that $\mathcal{U}\subset\mathbb{R}^{m}$ is a nonempty
compact convex set and contains the origin as its interior point.
$A\in\mathbb{R}^{n\times n}$ and $B\in\mathbb{R}^{n\times m}$ are unknown
system matrices with $(A,B)$ controllable, $(A,C)$ observable, satisfying
$|A|\leq A_{M}$, $|B|\leq B_{M}$.

Some standard assumptions are made on (1) and (2). Similar assumptions can be
found in [22, 26, 7] for solving (cooperative) output tracking problems.

Assumption 1. The input constraint has the box constrain form as
$\underline{u}\leq u\leq\overline{u}$ with elementwise inequality and
$\underline{u},\overline{u}$ the respective lower and upper bounds.

Assumption 2. There exists a constant matrix $K_{0}$ such that $A-BK_{0}$ is a
Hurwitz matrix with $-K_{0}x(t)\in\mathcal{U}$.

Assumption 3. The input relative degree (IRD) of system (1) is defined as
$\rho$.

Remark 1. In Assumption II, we refer the set $\mathcal{U}$ as a box constrain,
which accurately describes nearly any set of standard mechanical actuators.
Assumption II is made such that the initially feasibility can be achieved for
the system (1). Assumption II imposes IRD for the system (1), which is used
for simplify the solving of regulation problem [7].

For (1) and (2), the output of the system should track the given reference
$y_{\textrm{d}}(t),t\geq 0$. The tracking error can be given as
$e(t)=y(t)-y_{\textrm{d}}(t),t\geq 0$. To that end, a sampled-data MPC, which
is based on the repeated solution of an open-loop optimal control problem, is
provided in this paper. At each time instant $t=t_{k}$, the state $y(t_{k})$
is measured and then, the controller predicts the system behavior in the
future over a prediction horizon $T$ by minimizing a certain objective cost
function. The procedure is repeated at every sampling time instant $t_{k}$ for
$k=1,2,\ldots$

For the system (1) and (2), the cost function
$J(t_{k}):=J(x(t_{k}),y_{\textrm{d}}(\cdot),\hat{u}_{k}(\cdot))$, at time
$t_{k}$, is defined as

$\displaystyle J(x(t_{k}),y_{\textrm{d}}(s),\hat{u}_{k}(s))=$
$\displaystyle\int_{t_{k}}^{t_{k}+T}L(x(t_{k}),y_{\textrm{d}}(s),\hat{u}_{k}(s))\textrm{d}s$
$\displaystyle+F(y_{\textrm{d}}(t_{k}+T),\hat{y}(t_{k}+T))$ (3)

where, $L(\cdot,\cdot,\cdot)$ and $F(\cdot,\cdot)$ denote the stage and
terminal cost functions with the from

$\displaystyle
L(x(t_{k}),y_{\textrm{d}}(t),\hat{u}_{k}(t))=\|e(t)\|_{Q}^{2}+\|\hat{u}_{k}(t))\|_{R}^{2},\quad
t\geq 0,$

where $Q=Q^{T}\succ 0$ and $R=R^{T}\succeq 0$ are symmetric and sign definite
weight matrices.

Then, at time $t_{k}$, the optimal control signal ${\hat{u}}_{k}^{\star}(s)$,
$s\in[t_{k},t_{k}+T]$, is obtained by solving the following finite-horizon
optimal control problem as

$\displaystyle\min\limits_{\hat{u}_{k}(s)\in\mathcal{U}}J(x(t_{k}),y_{\textrm{d}}(s),\hat{u}_{k}(s))$
$\displaystyle
s.t.\quad{\dot{\hat{x}}}(s)=\mathcal{H}\left(\hat{x}(s),\hat{u}_{k}(s)\right){\Theta}$
(4a) $\displaystyle\qquad~{}\hat{y}(s)=C\hat{x}(s),{\hat{x}}(t_{k})=x(t_{k})$
(4b) $\displaystyle\qquad~{}{\hat{u}}(s)\in\mathcal{U},\quad
s\in[t_{k},t_{k}+T]$ (4c)

where
$\mathcal{H}(\cdot,\cdot):\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{n\times({n^{2}+mn})}$
is defined as

$\displaystyle\mathcal{H}(x,u)=$ $\displaystyle\left[(x\otimes
I_{n})^{\textrm{T}}~{}~{}(u\otimes I_{n})^{\textrm{T}}\right]$ (5)

and $\Theta$ denotes the vector of the dynamical parameters for (1), defined
as

$\displaystyle\Theta=\left[\text{vec}(A)^{\textrm{T}}~{}~{}\text{vec}(B)^{\textrm{T}}\right]^{\textrm{T}}\in\mathbb{R}^{n^{2}+nm}$
(6)

where $\text{vec}(\cdot)$ denotes the vectorization operator, defined as for
any $P=[p_{ij}]\in\mathbb{R}^{n\times l},l\in\\{n,m\\}$,
$\text{vec}(P)=\hat{P}$, where

$\displaystyle\hat{P}=[p_{11},p_{21},...,p_{n1},p_{12},...,p_{n2},...,p_{nl}]^{\textrm{T}}$
(7)

If the controller has the exact knowledge of system matrices $A$ and $B$, then
we can directly solve the finite-time optimal control problem (4) by [27] to
get the desired predictive control policy ${\hat{u}}_{k}^{\star}$. But, in
this paper, we seek a data-driven approach to remove the model-based
dependence of (6) on the predictive controller design. Thus, the primary
objective of this paper is to explore an online learning-based approach to
find the data-driven predictive control policy, without requiring any
knowledge of the system dynamics, just by using _a prior_ data of inputs and
outputs.

##  III Data-driven predictive control with completely unknown dynamics

In this section, to facilitate the predictive controller design, the state
${\hat{x}}$ and the system parameters $\Theta$ of the predicted model (4a) are
both estimated from the input-output measurements, using a simultaneous state
and parameter estimator. Our proposed online-learning strategy does not rely
on $\Theta$, i.e., either $A$ or $B$, which is totally data-driven approach.

###  III-A Dynamical Parameters Estimator

To obtain the dynamical parameters $\Theta$ in (6), i.e., $A$ and $B$, we
present our online learning strategy with the previous input-output data. By
using more prior data for continuous-time system (1) than the instants of
solving (4), we try to seek a least-squares optimization solution to get the
estimations of all the dynamical parameters $\Theta$. To that end, at each
time instant $t=t_{k}$, with enough prior date collected, we can find a
minimal periodic time interval $\delta t$, such that for some $t_{j}<t_{k}$ of
$j\in\\{1,2,...,k-1\\}$, there exists a constant integer $N_{k}$ satisfying
$N_{k}\delta t=t_{k}-t_{j}$.

To begin with, for any time instant $t\leq t_{k}$, doing the integrals at the
time interval $[t-\delta t,t]$ along the trajectory of (1) leads to

$\displaystyle x(t)-x(t-\delta t)=$ $\displaystyle A\int_{t-\delta
t}^{t}x\textrm{d}\tau+B\int_{t-\delta t}^{t}\hat{u}^{\star}\textrm{d}\tau$
$\displaystyle=$ $\displaystyle\mathcal{H}\left(\int_{t-\delta
t}^{t}x\textrm{d}\tau,\int_{t-\delta
t}^{t}\hat{u}^{\star}\textrm{d}\tau\right){\Theta}$ (8)

By rearranging (8), we have the linear error system in the form of

$\displaystyle\mathcal{F}(t)=\mathcal{G}(t)\Theta,~{}~{}~{}\forall
t\in\mathbb{R}_{[t_{j},t_{k}]}$ (9)

where the matrices $\mathcal{F}(\cdot):\mathbb{R}_{\geq
0}\rightarrow\mathbb{R}^{n}$ and $\mathcal{G}(\cdot):\mathbb{R}_{\geq
0}\rightarrow\mathbb{R}^{n\times(n^{2}+m^{2})}$ are defined,

$\displaystyle\mathcal{F}(t)=$ $\displaystyle\begin{cases}x(t)-x(t-\delta
t),\quad t\in[\delta t,\infty)\\\ 0,\qquad\qquad\qquad\quad t<\delta
t\end{cases}$ (10) $\displaystyle\mathcal{G}(t)=$
$\displaystyle\left[(\mathcal{A}(t)\otimes
I_{n})^{\textrm{T}}~{}~{}(\mathcal{B}(t)\otimes I_{n})^{\textrm{T}}\right]$
(11)

where $\otimes$ denotes the Kronecker product, and the vectors
$\mathcal{A}(\cdot):\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}^{n}$ and
$\mathcal{B}(\cdot):\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}^{n}$ are defined,

$\displaystyle\mathcal{A}(t)=$ $\displaystyle\begin{cases}\int_{t-\delta
t}^{t}x(\tau)\textrm{d}\tau,\quad t\in[\delta t,\infty)\\\
0,\qquad\qquad\quad\quad t<\delta t\end{cases}$ (12)
$\displaystyle\mathcal{B}(t)=$ $\displaystyle\begin{cases}\int_{t-\delta
t}^{t}\hat{u}^{\star}(\tau)\textrm{d}\tau,\quad t\in[\delta t,\infty)\\\
0,\qquad\qquad\quad\quad t<\delta t\end{cases}$ (13)

By (12)-(13), we note that for any time instant $t\leq t_{k}$, by using the
previous input-output measurements $x(\tau)$ and $u(\tau)$, $\tau\leq t$,
$\mathcal{F}(t)$ and $\mathcal{G}(t)$ in (10)-(11) are both available. Then,
by taking periodic sampling $\delta t$ for (1) at each time $t=t_{k}-i\delta
t$ with $i\in\\{0,1,2,...,N_{k}\\}$, such that $t_{j}\leq t\leq t_{k}$, then
we can get a set of composite data,

$\displaystyle\mathfrak{D}_{k}=\bigcup\\{\mathcal{F}_{i},~{}\mathcal{G}_{i}\\}_{i=0}^{N_{k}}$
(14)

which satisfies

$\displaystyle\mathcal{F}_{i}=\mathcal{G}_{i}\hat{\Theta}+\epsilon,~{}\forall
i\in\\{1,2,...,N_{k}\\}$ (15)

where $\mathcal{F}_{i}=\mathcal{F}(t_{k}-i\delta t)$,
$\mathcal{G}_{i}=\mathcal{G}(t_{k}-i\delta t)$,
$\hat{\Theta}=[\hat{A}^{T}~{}~{}\hat{B}^{T}]^{T}$ is an estimate of $\Theta$
in (6) and $\epsilon$ is the estimation error due to the data-driven
approximation of $\Theta$ by using $\mathfrak{D}_{k}$.

To bring $\epsilon$ to its minimum value, we have the following optimization
problem (OP):

$\displaystyle\min\limits_{\hat{\Theta}}~{}~{}\epsilon^{\textrm{T}}\epsilon$
(16) $\displaystyle s.t.\quad(15)~{}\text{and}~{}(9)-(13)$

Then, an adaptive least-square method is used here to give a solution to
problem (16). But, before that, the full rank of dataset $\mathfrak{D}_{k}$ is
defined as follows. To guarantee the existence of the solution of (16), we
also give the full rank condition for dataset $\mathfrak{D}_{k}$.

Definition 1.  At each time instant $t_{k}$ for some integer $N_{k}$, the data
stack $\mathfrak{D}_{k}$ is said to have full rank, if there exists an integer
$N_{0}>0$, such that for all $N_{k}\geq N_{0}$, we have the matrix
$\mathfrak{A}_{k}$, defined as
$\mathfrak{A}_{k}:=\sum_{i=1}^{N_{k}}\mathcal{G}_{i}^{\textrm{T}}\mathcal{G}_{i}\in\mathbb{R}^{(n^{2}+nm)\times(n^{2}+nm)}$,
satisfying

$\displaystyle\text{rank}(\mathfrak{A}_{k})=n^{2}+nm$ (17)

Then, we have the following lemma to give a sufficient condition to guarantee
the full rank property of dataset $\mathfrak{D}_{k}$. The proof is similar to
[28], and here it is omitted.

Lemma 1. A data stack $\mathfrak{D}_{k}$ has full rank, if there exists a
constant $\underline{d}>0$, such that

$\displaystyle 0<\underline{d}<\gamma_{m}(\mathfrak{A}_{k})$ (18)

Assume that the measurements of the inputs and the outputs are prior collected
at an enough large number $N_{k}\gg n^{2}+nm$ points of time $t_{k}-i\delta t$
with $i\in\\{0,1,...,N_{k}\\}$, which makes data stack $\mathfrak{D}_{k}$ has
full rank by Lemma III-A. Then, at time $t=t_{k}$, one can use the data stack
$\mathfrak{D}_{k}$ to evaluate the unknown dynamical parameters $\Theta$ in
(6). By solving (16), the following learning-based update law is obtained

$\displaystyle\dot{\hat{\Theta}}=\eta_{\theta}\sum\limits_{i=1}^{N_{k}}\mathcal{G}_{i}^{\textrm{T}}\left(\mathcal{F}_{i}-\mathcal{G}_{i}\hat{\Theta}\right)$
(19)

where $0<\eta_{\theta}\in\mathbb{R}$ is a constant learning rate.

For ease of exposition, the original continuous-time system (1) is expressed
in the form of,

$\displaystyle{\dot{x}}(t)=\mathcal{H}\left(x(t),u(t)\right){\Theta}$ (20)

By considering the linear error system with parameters update law (19), it
follows that

$\displaystyle{\dot{x}}(t)=$
$\displaystyle\mathcal{H}\left(x(t),u(t)\right)\hat{\Theta}+w(x(t),u(t))$
$\displaystyle:=$ $\displaystyle\tilde{A}x(t)+\tilde{B}u(t)+w(t)$ (21)

where $\tilde{B}=\text{vec}^{-1}(\hat{B})$,
$\tilde{A}=\text{vec}^{-1}(\hat{A})$, $\text{vec}^{-1}(\cdot)$ denotes the
converse vectorization operator, that is, for any vector
$\hat{P}\in\mathbb{R}^{n\times l}$ in (7), we have
$\text{vec}^{-1}(\hat{P})=\tilde{P}=[p_{ij}]\in\mathbb{R}^{n\times l}$. And
$w(t):=w(x(t),u(t))$ is the continuous approximation error resulting from
parameters uncertainty (19). If $\epsilon=0$, then $\hat{\Theta}=\Theta$, it
implies $\tilde{A}=A$ and $\tilde{B}=B$, thus we have $w(x,u)=0$.

The following theorem analyses the property of the term $w$ with respect to
(21).

Theorem 1. The approximate error $w(t)$ in (20) is slowly time-varying,
bounded and satisfies $\lim_{t\rightarrow\infty}w(t)=0$.
_Proof._ By considering the closed-loop dynamics (20) and (21), we can refer
to $w(t)$ as the unknown disturbance caused by the parameters uncertainty of
(19). Then, $w(t)$ satisfies

$\displaystyle
w(t)=\mathcal{H}\left(x(t),u(t)\right)\left(\hat{\Theta}-\Theta\right)$ (22)

Letting $\tilde{\Theta}=\hat{\Theta}-\Theta$ and bringing (19) to (22) lead to

$\displaystyle|\dot{w}(t)|\leq$
$\displaystyle|\mathcal{H}\left(x(t),u(t)\right)||\dot{\tilde{\Theta}}|$
$\displaystyle\leq$
$\displaystyle\eta_{\theta}|\mathcal{H}\left(x(t),u(t)\right)||\sum\limits_{i=1}^{N_{k}}\mathcal{G}_{i}^{\textrm{T}}\mathcal{G}_{i}||\dot{\tilde{\Theta}}|$
$\displaystyle=$
$\displaystyle\eta_{\theta}|\mathcal{H}\left(x(t),u(t)\right)||\mathfrak{A}_{k}||\dot{\tilde{\Theta}}|$
(23)

By Lemma 1, we have
$\gamma_{m}(\mathfrak{A}_{k})\leq|\mathfrak{A}_{k}|\leq\gamma_{M}(\mathfrak{A}_{k})$.
Then, by (19), it implies
$\dot{\tilde{\Theta}}=-\eta_{\theta}\mathfrak{A}_{k}\tilde{\Theta}$ and thus
$\lim_{t\rightarrow\infty}\tilde{\Theta}(t)=0$ and $\tilde{\Theta}\leq
e^{-\eta_{\theta}\gamma_{m}(\mathfrak{A}_{k})}\Theta(t_{k}-N_{k}\delta t)$.
So,
$|\dot{\tilde{\Theta}}|\leq\eta_{\theta}\gamma_{M}(\mathfrak{A}_{k})e^{-\eta_{\theta}\gamma_{m}(\mathfrak{A}_{k})}\Theta(t_{k}-N_{k}\delta
t):=C_{\Theta}$. Besides, for a fixed dataset $\mathfrak{D}_{k}$ in (14), we
have $|\mathcal{H}\left(x(t),u(t)\right)|$ in (23) bounded for some real
constant $M\in\mathbb{R}$ at each time $t=t_{k}-i\delta t$. Thus, we have

$\displaystyle|\dot{w}(t)|\leq\eta_{\theta}MC_{\Theta}\gamma_{M}(\mathfrak{A}_{k})$
(24)

It implies that $w(t)$ is slowly time-varying and
$\lim_{t\rightarrow\infty}w(t)=0$. By (23),(24), we also have $|w(t)|\leq
e^{\eta_{\theta}MC_{\Theta}\gamma_{M}(\mathfrak{A}_{k})}$, which means that
$w(t)$ is bounded.

Remark 2.  Note that in (III-A), the conventional receding horizon expression
$Ax(t)+Bx(t)$ depending on the unknown matrices $A,B$ is replaced by the term
$\mathcal{H}\left(\hat{x}(s),\hat{u}_{k}(s)\right)\tilde{\Theta}$, where
$\tilde{\Theta}$ can be obtained by repeatedly learning from the states and
inputs measurements. Furthermore, this learned results will not affects the
convergence of the system by Lemma III-A. Therefore, (19) plays an important
role in identifying the system dynamics from the _a prior_ data. As a result,
the requirement of the system matrices in predicting the behavior of (II) can
be replaced by the state and input information measured online.

###  III-B Receding-horizon Optimization

To facilitate the data-driven predictive controller design for the system (21)
with the dynamical parameters estimator (19), the receding-horizon predictive
control problem of (4), at time instant $t=t_{k}$, can be reformulated,

$\displaystyle{\hat{u}}_{k}^{\star}(s)=\arg\mathop{\min}\limits_{{{\hat{u}}}(t)\in\mathcal{U}}{J}({{x}}({t_{k}}),y_{\textrm{d}}(s),{{\hat{u}_{k}}}(s))$
$\displaystyle
s.t.\quad{\dot{\hat{x}}}(s)=\tilde{A}{\hat{x}}(s)+\tilde{B}\hat{u}_{k}(s),$
(25a)
$\displaystyle\qquad~{}\hat{y}(s)=C\hat{x}(s),{\hat{x}}(t_{k})=x(t_{k}),$
(25b) $\displaystyle\qquad~{}{\hat{u}}(s)\in\mathcal{U},\quad
s\in[t_{k},t_{k}+T].$ (25c)

To solve the optimization problem (25), under Assumption II, we define the
decision variables as
${\bar{u}}_{k}(s)=[\hat{u}_{k}^{\textrm{T}}(s),(\hat{u}_{k}^{[1]})^{\textrm{T}}(s),\ldots,(\hat{u}_{k}^{[r]})^{\textrm{T}}(s)]$
for some control order $r$ larger than $\rho\geq 1$. Note that the first term
of ${\bar{u}}_{k}(\tau)$ is the to-be-optimized control input
$\hat{u}_{k}(\tau)$ in (25). More generally, for the control law
$\hat{u}_{k}(\tau)$ with a large enough control order $r$, we let
$\hat{u}_{k}^{[l]}(\tau)=0$ for any integer $l\geq r$.

Then, for the output prediction of optimization problem (25), by following
(2), the future output $y(s)=y(t+\tau),t=t_{k},k=1,2,\ldots$, in the moving
horizon $\tau\in[0,T]$ is approximated by Taylor series expansion,

$\displaystyle y(t+\tau)=y(t)+\tau
y^{[1]}(t)+\cdots+\frac{\tau^{r}}{r!}y^{[r]}(t)+O(\tau^{r})$ (26)

where the $i$-th derivative of the output $y^{[i]}(t)$ with
$i\in\\{1,2,...,\rho,...,r\\}$ is obtained by

$\displaystyle{y^{[i]}}=$ $\displaystyle
C{{\tilde{A}}^{i}}x+\sum\limits_{k=0}^{i-1}C{{\tilde{A}}^{i-1-k}}{w^{[k]}},i=1,...,\rho-1$
(27) $\displaystyle{y^{[j]}}=$ $\displaystyle
C{{\tilde{A}}^{j}}x+\sum\limits_{k=0}^{j-\rho}C{{\tilde{A}}^{j-1-k}}\tilde{B}{u^{[k]}}$
$\displaystyle+\sum\limits_{k=0}^{j-1}C{{\tilde{A}}^{j-1-k}}{w^{[k]}},j=\rho,...,r$
(28)

By rewriting the output $y(t+\tau)$ in a compact form, it follows that

$\displaystyle
y(t+\tau)=\begin{bmatrix}T_{1}(\tau)&T_{2}(\tau)\end{bmatrix}\begin{bmatrix}Y_{1}\\\
Y_{2}\end{bmatrix}$ (29)

where
${{T}}_{1}(\tau)=\left[1,\tau,\ldots,\frac{\tau^{\rho-1}}{(\rho-1)!}\right]$,
$T_{2}(\tau)=\left[\frac{\tau^{\rho}}{\rho!},\ldots,\frac{\tau^{r}}{r!}\right]$,
${Y}_{1}=\left[y^{\textrm{T}},(y^{[1]})^{\textrm{T}},\ldots,(y^{[\rho-1]})^{\textrm{T}}\right]^{\textrm{T}}$,
${Y}_{2}=\left[(y^{[\rho]})^{\textrm{T}},(y^{[\rho+1]})^{\textrm{T}},\ldots,(y^{[r]})^{\textrm{T}}\right]^{\textrm{T}}$,
and

$\displaystyle Y_{1}=$ $\displaystyle\mathcal{A}_{1}x+\mathcal{B}_{1}\bar{w}$
(30) $\displaystyle Y_{2}=$
$\displaystyle\mathcal{A}_{2}x+\mathcal{B}_{2}\bar{w}+\mathcal{B}_{3}\bar{u}$
(31)

where
${\bar{w}}(s)=\left[w^{\textrm{T}}(s),(w^{[1]})^{\textrm{T}}(s),\ldots,(w^{[r]})^{\textrm{T}}(s)\right]^{\textrm{T}}$
and $\mathcal{A}_{1},\mathcal{B}_{1},\mathcal{A}_{2},\mathcal{B}_{2}$ and
$\mathcal{B}_{3}$ are defined,

$\displaystyle\mathcal{A}_{1}=$ $\displaystyle\begin{bmatrix}C\\\
C\tilde{A}\\\ \vdots\\\
C\tilde{A}^{\rho-1}\end{bmatrix},\mathcal{B}_{1}=\begin{bmatrix}C&0&\cdots&0\\\
C\tilde{A}&C&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
C\tilde{A}^{\rho-1}&C\tilde{A}^{\rho-2}\tilde{B}&\cdots&C\end{bmatrix}$
$\displaystyle\mathcal{A}_{2}=$
$\displaystyle\begin{bmatrix}C\tilde{A}^{\rho}\\\ C\tilde{A}^{\rho+1}\\\
\vdots\\\
C\tilde{A}^{r}\end{bmatrix},\mathcal{B}_{2}\\!=\\!\begin{bmatrix}C\tilde{A}^{\rho-1}&0\cdots&0\\\
C\tilde{A}^{\rho}&C\tilde{A}^{\rho-1}&\cdots&0\\\
\vdots&\vdots&\ddots&\cdots\\\
C\tilde{A}^{r}&C\tilde{A}^{r-1}&\cdots&C\tilde{A}^{\rho-1}\end{bmatrix}$
$\displaystyle\mathcal{B}_{3}=$
$\displaystyle\begin{bmatrix}C\tilde{A}^{\rho-1}\tilde{B}&0&\cdots&0\\\
C\tilde{A}^{\rho}\tilde{B}&C\tilde{A}^{\rho-1}\tilde{B}&\cdots&0\\\
\vdots&\vdots&\ddots&\cdots\\\
C{\tilde{A}}^{r}\tilde{B}&C{\tilde{A}}^{r-1}\tilde{B}&\cdots&C{\tilde{A}}^{\rho-1}\tilde{B}\end{bmatrix}$
(32)

Then, for the reference signal $y_{d}(t+\tau)$, by Theorem 1, we have
$w^{[k]}(\tau)\approx 0$ for $k=1,2,\ldots,r$. Thus, $y_{d}(t+\tau)$ satisfies
(29)-(31) with
$Y_{1,d}=\left[y_{d}^{\textrm{T}},(y_{d}^{[1]})^{\textrm{T}},\ldots,(y_{d}^{[\rho-1]})^{\textrm{T}}\right]$
and
$Y_{2,d}=\left[(y_{d}^{[\rho]})^{\textrm{T}},(y_{d}^{[\rho+1]})^{\textrm{T}},\ldots,(y_{d}^{[r]})^{\textrm{T}}\right]$.
But, along with the solution of optimization problem (25) without considering
the parameters uncertainty from (21), it leads to

$\displaystyle\hat{y}(t+\tau)=\begin{bmatrix}T_{1}(\tau)&T_{2}(\tau)\end{bmatrix}\begin{bmatrix}Y_{1}\\\
Y_{2}\end{bmatrix}$ (33) $\displaystyle Y_{1}=\mathcal{A}_{1}x,\quad
Y_{2}=\mathcal{A}_{2}x+\mathcal{B}_{3}{\bar{u}}$ (34)
$\displaystyle\hat{u}(t+\tau)=T(\tau)\bar{u}$ (35)

where $T(\tau)=[T_{1}(\tau)~{}T_{2}(\tau)]$.

Then, by bring (34)-(35) into the given performance index (3), one can go to,

$\displaystyle J(t_{k})=$
$\displaystyle\int_{0}^{T}\left[\tilde{Y}_{1}^{\textrm{T}},\tilde{Y}_{2}^{\textrm{T}}\right]\begin{bmatrix}\Xi_{1}^{\textrm{T}}(\tau)\\\
\Xi_{2}^{\textrm{T}}(\tau)\end{bmatrix}\left[\Xi_{1}(\tau),\Xi_{2}(\tau)\right]\begin{bmatrix}\tilde{Y}_{1}\\\
\tilde{Y}_{2}\end{bmatrix}$
$\displaystyle+\int_{0}^{T}\hat{u}^{\textrm{T}}T_{3}^{\textrm{T}}(\tau)RT_{3}(\tau)\bar{u}\textrm{d}\tau+F(\tilde{Y}_{i}(t_{k}+T))$
(36)

where $\tilde{Y}_{i}=Y_{i}-Y_{i,d}$ for $i\in\\{1,2\\}$, and
$\Xi_{i}(\tau)=\sqrt{Q}T_{i}(\tau)$ for $i\in\\{1,2\\}$. By defining
${\mathcal{T}}_{i,j}=\int_{0}^{T}\Xi_{i}^{\textrm{T}}\Xi_{j}\textrm{d}\tau$
with $i,j\in\\{1,2\\}$ and
${\mathcal{T}}=\int_{0}^{T}T^{\textrm{T}}RT\textrm{d}\tau$, we note that
${\mathcal{T}}_{1,2}={\mathcal{T}}_{2,1}$. Thus, the performance index (37)
can also be rewritten as

$\displaystyle J(t_{k})=$
$\displaystyle\tilde{Y}_{1}^{\textrm{T}}{\mathcal{T}}_{1,1}\tilde{Y}_{1}+2\tilde{Y}_{1}^{\textrm{T}}{\mathcal{T}}_{1,2}\tilde{Y}_{2}+\tilde{Y}_{2}^{\textrm{T}}{\mathcal{T}}_{2,2}\tilde{Y}_{2}$
$\displaystyle+{\bar{u}}^{\textrm{T}}\mathcal{T}{\bar{u}}+F(\tilde{Y}_{i}(t_{k}+T))$
(37)

Due to equation (34), taking partial derivative of $J(t_{k})$ with respect to
${\bar{u}}_{k}$ yields

$\displaystyle\frac{{\partial J}}{{\partial\bar{u}}}=$ $\displaystyle
2\left((\frac{{\partial\tilde{Y}_{2}}}{{\partial\bar{u}}})^{\textrm{T}}{\mathcal{T}}_{2,2}\tilde{Y}_{2}+{\mathcal{T}}_{4}\right){\bar{u}}+2(\frac{{\partial\tilde{Y}_{2}}}{{\partial\bar{u}}})^{\textrm{T}}{\mathcal{T}}_{1,2}^{\textrm{T}}\tilde{Y}_{1}$
$\displaystyle=$ $\displaystyle
2\left({\mathcal{B}}_{3}^{\textrm{T}}{\mathcal{T}}_{2,2}{\mathcal{B}}_{3}+{\mathcal{T}}_{4}\right){\bar{u}}+2{\mathcal{B}}_{3}^{\textrm{T}}{\mathcal{T}}_{1,2}^{\textrm{T}}\tilde{Y}_{1}$
$\displaystyle+2{\mathcal{B}}_{3}^{\textrm{T}}{\mathcal{T}}_{2,2}\left(\mathcal{A}_{2}x-Y_{2,d}\right)$
(38)

By letting ${\partial J}/{{\partial{\bar{u}}_{k}}}=0$, we can get the
optimized predictive control law ${\bar{u}}_{k}^{\star}$ as

$\displaystyle{\bar{u}}_{k}^{\star}=$
$\displaystyle-\left({\mathcal{B}}_{3}^{\textrm{T}}{\mathcal{T}}_{2,2}{\mathcal{B}}_{3}+{\mathcal{T}}_{4}\right)^{-1}{\mathcal{B}}_{3}^{\textrm{T}}\left({\mathcal{T}}_{2,2}\bar{Y}_{2}-{\mathcal{T}}_{1,2}^{\textrm{T}}\tilde{Y}_{1}\right)$
(39)

where $\bar{Y}_{2}=Y_{2,d}-\mathcal{A}_{2}x$. Then, taking the first row of
the optimized control law (39), the continuous-time predictive control law,
applied to the plant, is given by

$\displaystyle\hat{u}_{k}^{\star}(t)=I_{u}{\bar{u}}_{k}^{\star}$ (40)

where $I_{u}=[1,0,...,0]_{1\times(r+1)}$.

Remark 3. Note that from (40), the existence of the optimized solution
${\bar{u}}_{k}^{\star}$ depends on the reversibility of matrix
$\mathcal{M}={\mathcal{B}}_{3}^{\textrm{T}}{\mathcal{T}}_{2,2}{\mathcal{B}}_{3}+{\mathcal{T}}_{4}$,
with ${\mathcal{B}}_{3}$ computed from (32), $\tilde{A}$ and $\tilde{B}$
calculated from (19). Thus, before we implement the receding-horizon
optimization, we first check the reversibility of matrix $\mathcal{M}$ by
removing the repeated columns of dataset $\mathfrak{D}_{k}$, only left the
distinct columns.

###  III-C Handing constraints

To deal with optimal control problems, MPC can allow for industrial processes
uncertainties and constraints much more straightforwardly than other methods
[2], [3]. Assumption II will allow us to use a specialized active-set
algorithm which is more efficient and easier to implement. A box-constraint
solver can be immediately generalized to any linear inequality constraints
using slack variables. In the following, we formalize two classical ways to
enforce the control limits.

####  III-C1 Saturating Functions

A conventional attempt to enforce box constraints is to clamp the controls in
the forward-pass. The element-wise clamping, or projection operator, is
denoted by $sat\text{($\cdot$)}$,which is the input saturation function
defined as

$\displaystyle sat(u)$
$\displaystyle=\left[sat(u_{1})\>sat(u_{2})\>\cdots\>sat(u_{m})\right]$ (41)
$\displaystyle sat(u_{i}(t))$ $\displaystyle=\begin{cases}u_{i}(t)&{\textstyle
if}\,u_{i,min}<u_{i}(t)<u_{i,max}\\\ u_{i,min}&{\displaystyle
if}\,u_{i}(t)\leq u_{i,min}\\\ u_{i,max}&if\,u_{i}(t)\geq
u_{i,max}\end{cases}$

with $u=[u_{1}\>u_{2}\cdots u_{m}]$, and $u_{i,min}\leq 0$ and $u_{i,max}\geq
0$ are the boundaries of $i$th control input of system (1).

It implies that it is tempting to simply replace the obtained control in the
forward-pass with

$\hat{u}=sat(u^{*})$

However, the corresponding search direction may not be a descent direction
anymore, harming convergence.

####  III-C2 Squashing Functions

Another way to enforce box constraints is to introduce a sigmoidal squashing
function $s(u)$ on the controls

$x_{i+1}=f(x_{i},s(u_{i}))$

where $s(\cdot)$ is an element-wise sigmoid with the vector limits

$\underset{u\to-\infty}{lim}s(u)=\underline{u},\quad\underset{u\to\infty}{lim}s(u)=\overline{u}$

For example,
$s(u)=\frac{\overline{u}-\underline{u}}{2}tanh(u)+\frac{\overline{u}+\underline{u}}{2}$
is such a function. A cost term should be kept on the original $u$ and not
only on the squashed $s(u)$. Otherwise it will reach very high or low values
and get stuck on the plateau. An intuition for the poor practical performance
of squashing is given by the nonlinearity of the sigmoid. Since the backward
pass uses a locally quadratic approximation of the dynamics, significant
higher order terms will always have a detrimental effect on convergence.

###  III-D Our Proposed Algorithm

Our proposed method is summarized as the following Algorithm $1$.

Data: current period $t$; a initially stable control $u^{\star}=-K_{0}x$;
prediction horizon $T$; terminal cost $F$;

Result: Optimal MPC input $u^{\ast}_{t}$

1 $t_{k}=k\leftarrow 0$;

2 Collect data $\mathfrak{D}_{k}$ in (14) ;

3 while  __ do

4 Collect the data and form $\mathfrak{D}_{k}$;

5 if  _(18) is satisfied_ then

6 Generate the estimator (19) by using $\mathfrak{D}_{k}$;

7 Implement (25) to obtain the optimized control (40);

8 Time evolves continuously with $t$;

9 end if

10 $k\leftarrow k+1$;

11 end while

Algorithm 1 Data-driven predictive control algorithm

###  III-E Performance Analysis

Before proceeding further, we first introduce the following definition and
lemma.

Definition 2. For the system (21), given a compact set $\mathbb{E}$, with
$\\{0\\}\subset\mathbb{E}\subseteq\mathbb{R}^{n}$ and $\mathbb{E}$ being a
robustly positively invariant set, if there exists a positive definite
function $V(\cdot):\mathbb{R}^{n}\to\mathbb{R}_{\geq 0}$, such that,

$\displaystyle V(x)$
$\displaystyle\geq\alpha_{1}(|x|),~{}~{}V(x)\leq\alpha_{2}(|x|)+c_{1}$ (42)
$\displaystyle V(\dot{x})$
$\displaystyle-V(x)\leq-\alpha_{3}(|x|)+\alpha_{4}(|w|)+c_{2}$ (43)

for all $t\in\mathbb{R}_{\geq 0}$ with $\alpha_{1},\alpha_{2},\alpha_{3}$
being $\mathcal{K}_{\infty}$ function, $\alpha_{4}$ being $\mathcal{K}$
function, and $c_{1},c_{2}\geq 0$. Then, the function $V(\cdot)$ is a regional
input-to-state practical stability (ISpS)-type Lyapunov function in
$\mathbb{E}$ for the system.

Based on Definition 2, we can have the following lemma, directly borrowed from
[7].

Lemma 2. Given a robust positively invariant set $\mathbb{E}$ for the system
(21), if it admits an ISpS-type Lyapunov function $V(\cdot)$, then the system
is regional ISpS in $\mathbb{E}$, and all the signals of the closed-loop
system with the control input $u_{\mathbb{E}}$ are bounded, where
$u_{\mathbb{E}}(t)$ denotes the control such that the set $\mathbb{E}$ is an
invariant region satisfying the constraints.

Note that by Assumption 1, the developed optimization problem of data-driven
predictive control in (25) is initially feasible with $u_{0}=-K_{0}x$, then
the global stability can be proved by using the Lemma 2.

Theorem 2. Suppose the Assumption 1 and 2 hold for system with a robust
positively invariant set $\mathbb{E}$, then the closed-loop control plant (1)
under the continuous-time MPC law (40), is globally asymptotically stable.

_Proof._ The proof is composed by two parts, feasibility and convergence.

_Feasibility:_ Consider any time $t_{k}$ such that the problem of (25) has a
solution and the optimal input $\hat{u}_{k}^{\star}$ is implemented for time
$[t_{k},t_{k+1})$. Assumed that at $t_{k+1}$, $\hat{y}(t_{k+1})=y(t_{k+1})$.
Therefore, the remaining piece of optimal input
$\hat{u}_{k}^{\star}(s),\tau\in[t_{k+1},t_{k}+T]$ satisfies the input
constraints. Thus, we construct the control input as,

$\displaystyle\hat{u}_{k+1}(\tau)=\begin{cases}\hat{u}_{k}^{\star}(\tau),~{}~{}\tau\in[t_{k+1},t_{k}+T]\\\
u_{\mathbb{E}}(\tau),~{}~{}\tau\in[t_{k}+T,t_{k+1}+T]\end{cases}$ (44)

where $u_{\mathbb{E}}(\tau)$ makes the desired reference reached and the
constraints satisfied. Thus, the predictive control problem is feasible at
$t_{k+1}$. It implies that feasibility of the problem at $t_{k}$ implies the
recursive feasibility at $t_{k+1}$.

_Convergence:_ Let the optimal cost function at $t_{k}$ as the value function
$V(x(t_{k}))=J^{\star}(x(t_{k}),y_{\textrm{d}}(s),\hat{u}_{k}^{\star}(s))$. If
$V(x(t_{k}))$ is strictly decreasing, the tracking error $e$ will converge to
the origin. To this end, we write the value function at $t_{k}$ as,

$\displaystyle V(t_{k})=$
$\displaystyle\int_{t_{k}}^{t_{k}+T}\left(\|e(\tau)\|_{Q}^{2}+\|\hat{u}_{k}^{\star}(\tau))\|_{R}^{2}\right)\textrm{d}\tau$
$\displaystyle+F(y_{d}(t_{k}+T),\hat{y}(t_{k}+T))$ (45)

Then, by applying $\hat{u}_{k+1}(t)$ in (44) to the system, beginning from
$y(t_{k+1})$, one has

$\displaystyle J(t_{k+1})=$
$\displaystyle\int_{t_{k+1}}^{t_{k+1}+T}\left(\|e(\tau)\|_{Q}^{2}+\|\hat{u}_{k+1}(\tau))\|_{R}^{2}\right)\textrm{d}\tau$
$\displaystyle+F(y_{\textrm{d}}(t_{k+1}+T),\hat{y}(t_{k+1}+T))$ (46)

By substituting (45) in (46), it follows,

$\displaystyle J(t_{k+1})=$ $\displaystyle
V(t_{k})-\int_{t_{k}}^{t_{k+1}}\left(\|e(\tau)\|_{Q}^{2}+\|\hat{u}_{k}^{\star}(\tau))\|_{R}^{2}\right)\textrm{d}\tau$
$\displaystyle-F(y_{\textrm{d}}(t_{k}+T),\hat{y}(t_{k}+T))$
$\displaystyle+\int_{t_{k}+T}^{t_{k+1}+T}\left(\|e(\tau)\|_{Q}^{2}+\|\hat{u}_{k+1}(\tau))\|_{R}^{2}\right)\textrm{d}\tau$
$\displaystyle+F(y_{\textrm{d}}(t_{k+1}+T),\hat{y}(t_{k+1}+T))$ (47)

Note that $\mathbb{E}$ is a robust positively invariant set, it implies that
for all $x\in\mathbb{E}$, we have

$\displaystyle(\partial F/\partial
x)(Ax+Bu_{\mathbb{E}})+L(x,y_{\textrm{d}},u_{\mathbb{E}})\leq 0$ (48)

where
$L(x,y_{\textrm{d}},u_{\mathbb{E}})=\|e\|_{Q}^{2}+\|u_{\mathbb{E}}\|_{R}^{2}$.
By integrating the inequality (48) along the trajectory of
$\dot{\hat{x}}=\tilde{A}\hat{x}+\tilde{B}u_{\mathbb{E}}$ and
$\hat{y}=C\hat{x}$, we have
$\int_{t_{k}+T}^{t_{k+1}+T}\left(\|e(\tau)\|_{Q}^{2}+\|u_{\mathbb{E}}(\tau))\|_{R}^{2}\right)\textrm{d}\tau\leq
F(y_{\textrm{d}}(t_{k}+T),\hat{y}(t_{k}+T))-F(y_{\textrm{d}}(t_{k+1}+T),\hat{y}(t_{k+1}+T))$.
Thus, by (47), one gets,

$\displaystyle
V(t_{k})-J(t_{k+1})\\!\leq\\!-\int_{t_{k}}^{t_{k+1}}\left(\|e(\tau)\|_{Q}^{2}+\|\hat{u}_{k}^{\star}(\tau))\|_{R}^{2}\right)\textrm{d}\tau$
(49)

Further, based on $V(t_{k+1})=J^{\star}(t_{k+1})$, we obtain,

$\displaystyle
V(t_{k})-V(t_{k+1})\\!\leq\\!-\int_{t_{k}}^{t_{k+1}}\left(\|e(\tau)\|_{Q}^{2}+\|\hat{u}_{k}^{\star}(\tau))\|_{R}^{2}\right)\textrm{d}\tau$
(50)

It implies that $V(t_{k+1})$ is strictly decreasing. Hence, the proof is
complete.

$\displaystyle A=\begin{bmatrix}-17.98&-295.866&0&0&0&0\\\
0.0207&0.1889&0.0704&0&0&0\\\ 0&0.3879&0.8000&0&0&0\\\
0.0977&0&0&-18.01&-295.87&0\\\ 0&0.0617&0&0.0131&0.0433&0.0589\\\
0&0&0&0&0.3787&-0.622\end{bmatrix},\quad B=\begin{bmatrix}17.8996&-13.781\\\
-0.0131&0.0101\\\ 0&0\\\ 17.8636&17.8636\\\ 0.0082&0.0082\\\ 00\end{bmatrix}$
(50)

##  IV Application to two-CSTR process

Consider two continuous stirred tank reactor (CSTR) system with a full
description in [29, 30]. The open-loop model is a six-state continuous model.
The system matrices $A$ and $B$ are directly taken from [29], described in the
form of (1), as you can see in (50). The system output variables
$y_{1}=362.995x_{2}$ and $y_{2}=362.995x_{4}$, denoting the two tank outlet
temperatures. The control problem is to maintain the two tank temperatures at
desired values $y_{d}(t)=[y_{1d}(t)~{}y_{2d}(t)]^{\textrm{T}}$, where
$y_{1d}(t)=10$ when $0\leq t<5$s and $y_{1d}(t)=7$ when $t\geq 5$s,
$y_{2d}(t)=10$ when $0\leq t<5$s and $y_{2d}(t)=4$ when $t\geq 5$s. The
constraints is,

$\displaystyle\mathcal{U}=\\{u=[u_{1}~{}u_{2}]^{\textrm{T}}:|u_{1}|\leq
80,~{}~{}|u_{2}|\leq 70\\}$

In order to illustrate the efficiency of the proposed approach, the precise
knowledge of $A$ and $B$ is not used in the design of the predictive
controllers. Since the physical system is not stable, the initial stabilizing
feedback gain is set as $K_{0}$,

$\displaystyle
K_{0}=\begin{bmatrix}\begin{smallmatrix}-4.8949&-3426.8&-158.1712&-0.0320&-43.7963&-1.4675\\\
0.1&0&86.2934&1.1730&2.3886&104.8756\end{smallmatrix}\end{bmatrix}$

The weighting matrices $Q$ and $R$ are set to be
$Q=diag([10~{}100~{}10~{}10~{}100~{}10])$ and $R=diag([1~{}1])$, respectively.
In the simulation, the initial values are selected at the origin. The state
and input information is collected over each interval of $0.01$s. When time
arrives at $t=2$ s, all the inputs and outputs are repeatedly used to
approximate the matrices $A$ and $B$ with $\eta_{\theta}=0.85$. The predictive
control also starts at $t=2$ s with the prediction horizon $T=1$ s. Since
then, the control input is immediately updated by solving the problem (25),
and the convergence of of $A_{k}:=\tilde{A}$ and $B_{k}:=\tilde{B}$ to their
actual values is attained after 10 iterations. The procedure of solving (25)
is repeated over a fixed interval of $0.1$s. The convergence of $A_{k}$ and
$B_{k}$ to their actual values is illustrated in Fig.1. The trajectories of
the output variables and the flow rates are shown in Fig.2. It can be seen
that the data-driven predictive control algorithm can stabilize the system,
without requiring the system matrices.

Figure 1: Convergence of $A_{k}$ and $B_{k}$ to their actual values during the
control process.
Figure 2: The trajectories of the output variables and the flow rates.

##  V Conclusion

In this paper, a data-driven predictive control approach for continuous-time
linear system with completely unknown dynamics has been provided. This method
solves the infinite-horizon optimal control problem, using the system inputs
and outputs information collected online, without knowing the system matrices.
The methodology developed in this paper may serve as a computational tool to
study the finite-horizon adaptive optimal control of uncertain nonlinear
systems. Some related work has appeared in [15], which was developed using
neural networks, and also in our recent work [26], which proposed a framework
of distributed MPC to handle the asynchronous communication, using the _a
prior_ information associated with the interconnected neighbors to a
distributed optimal design.

## References

  * [1] P. Karamanakos, T. Geyer, N. Oikonomou, F. D. Kieferndorf, and S. Manias, “Direct model predictive control: A review of strategies that achieve long prediction intervals for power electronics,” _IEEE Industrial Electronics Magazine_ , vol. 8, no. 1, pp. 32–43, 2014.
  * [2] M. V. Kothare, V. Balakrishnan, and M. Morari, “Robust constrained model predictive control using linear matrix inequalities,” _Automatica_ , vol. 32, no. 10, pp. 1361–1379, 1996.
  * [3] D. Q. Mayne, J. B. Rawlings, C. V. Rao, and P. O. Scokaert, “Constrained model predictive control: Stability and optimality,” _Automatica_ , vol. 36, no. 6, pp. 789–814, 2000.
  * [4] Y. Xi, D. Li, and S. Lin, “Model predictive control - status and challenges,” _Acta Automatica Sinica_ , vol. 39, no. 3, pp. 222–236, 2013\.
  * [5] S. J. Qin, “Survey on data-driven industrial process monitoring and diagnosis,” _Annual Reviews in Control_ , vol. 36, no. 2, pp. 220–234, 2012\.
  * [6] Y. Xu, D. Li, Y. Xi, J. Lan, and T. Jiang, “Improved predictive controller on FPGA by hardware matrix inversion,” _IEEE Transactions on Industrial Electronics_ , 2018.
  * [7] J. Yang, W. X. Zheng, S. Li, B. Wu, and M. Cheng, “Design of a prediction-accuracy-enhanced continuous-time MPC for disturbed systems via a disturbance observer,” _IEEE Transactions on Industrial Electronics_ , vol. 62, no. 9, pp. 5807–5816, 2015.
  * [8] H. Shah and M. Gopal, “Model-free predictive control of nonlinear processes based on reinforcement learning,” in _IFAC-PapersOnLine_ , vol. 49, no. 1, 2016, pp. 89–94.
  * [9] X. Wang, B. Huang, and T. Chen, “Data-driven predictive control for solid oxide fuel cells,” _Journal of Process Control_ , vol. 17, no. 2, pp. 103–114, 2007.
  * [10] Z. Hou and S. Jin, “Data-driven model-free adaptive control for a class of MIMO nonlinear discrete-time systems,” _IEEE Transactions on Neural Networks_ , vol. 22, no. 12, pp. 2173–2188, 2011.
  * [11] D. Laurí, J. A. Rossiter, J. Sanchis, and M. Martínez, “Data-driven latent-variable model-based predictive control for continuous processes,” _Journal of Process Control_ , vol. 20, no. 10, pp. 1207–1219, 2010.
  * [12] D. Zhou, D. Sun, Z. Pang, and G. Liu, “Data-based predictive control for networked non-linear systems with two-channel packet dropouts,” _IET Control Theory & Applications_, vol. 9, no. 7, pp. 1154–1161, 2015.
  * [13] W. Aangenent, D. Kostic, B. de Jager, R. van de Molengraft, and M. Steinbuch, “Data-based optimal control,” in _Proceedings of the American Control Conference_ , 2005, pp. 1460–1465.
  * [14] S. S. Ge, Z. Li, and H. Yang, “Data driven adaptive predictive control for holonomic constrained under-actuated biped robots,” _IEEE Transactions on Control Systems Technology_ , vol. 20, no. 3, pp. 787–795, 2012.
  * [15] S. S. Ge, C. Yang, S. Member, and T. H. Lee, “Adaptive predictive control using neural network for a class of pure-feedback systems in discrete time,” _IEEE Transactions on Neural Networks_ , vol. 19, no. 9, pp. 1599–1614, 2008.
  * [16] R. Kadali, B. Huang, and A. Rossiter, “A data driven subspace approach to predictive controller design,” _Control Engineering Practice_ , vol. 11, no. 3, pp. 261–278, 2003.
  * [17] Y. Zhou, D. Li, Y. Xi, and Z. Gan, “Periodic event-triggered control for distributed networked multiagents with asynchronous communication: A predictive control approach,” _International Journal of Robust and Nonlinear Control_ , vol. 29, no. 1, pp. 43–66, 2019.
  * [18] G. Shi and R. E. Skelton, “Markov data-based lqg control,” _Journal of dynamic systems, measurement, and control_ , vol. 122, no. 3, pp. 551–559, 2000\.
  * [19] D. Li, Y. Xi, and F. Gao, “Synthesis of dynamic output feedback RMPC with saturated inputs,” _Automatica_ , vol. 49, no. 4, pp. 949–954, 2013.
  * [20] Y. Zhou, D. Li, Y. Xi, and Z. Gan, “Synthesis of model predictive control based on data-driven learning,” _Science China Information Sciences_ , vol. 63, pp. 1–3, 2020.
  * [21] D. Görges, “Relations between model predictive control and reinforcement learning,” in _IFAC-PapersOnLine_ , vol. 50, no. 1, 2017, pp. 4920–4928.
  * [22] Y. Jiang and Z. P. Jiang, “Computational adaptive optimal control for continuous-time linear systems with completely unknown dynamics,” _Automatica_ , vol. 48, no. 10, pp. 2699–2704, 2012.
  * [23] ——, _Robust Adaptive Dynamic Programming_. USA: Wiley-IEEE Press, 2017.
  * [24] Y. Zhou, K. G. Vamvoudakis, W. M. Haddad, and Z.-P. Jiang, “A secure control learning framework for cyber-physical systems under sensor and actuator attacks,” _IEEE Transactions on Cybernetics_ , 2020.
  * [25] ——, “A secure control learning framework for cyber-physical systems under sensor attacks,” in _2019 American Control Conference (ACC)_. IEEE, 2019, pp. 4280–4285.
  * [26] Y. Zhou, D. Li, J. Lu, Y. Xi, and L. Cen, “Networked and distributed predictive control of non-linear systems subject to asynchronous communication,” _IET Control Theory & Applications_, vol. 12, no. 4, pp. 504–514, 2018.
  * [27] F. L. Lewis, D. Vrabie, and V. L. Syrmos, _Optimal control_. John Wiley & Sons, 2012.
  * [28] P. Ioannou and J. Sun, _Robust adaptive control_. Upper Saddle River, NJ: PTR Prentice-Hall, 1996.
  * [29] Y. Cao and Z. Yang, “Multiobjective process controllability analysis,” _Computer and Chemical Engineering_ , vol. 28, pp. 83–90, 2004.
  * [30] T. Wang, H. Gao, and J. Qiu, “A combined fault-tolerant and predictive control for network-based industrial processes,” _IEEE Transactions on Industrial Electronics_ , vol. 63, no. 4, pp. 2529–2536, 2016.