# Quantum simulation of partial differential equations via Schrödingerisation

Shi Jin  Institute of Natural Sciences, School of Mathematical Sciences, MOE-
LSC, Shanghai Jiao Tong University, Shanghai, 200240, P. R. China  Shanghai
Artificial Intelligence Laboratory, Shanghai, China  Nana Liu
<EMAIL_ADDRESS>Institute of Natural Sciences, School of Mathematical
Sciences, MOE-LSC, Shanghai Jiao Tong University, Shanghai, 200240, P. R.
China  Shanghai Artificial Intelligence Laboratory, Shanghai, China
University of Michigan-Shanghai Jiao Tong University Joint Institute, Shanghai
200240, China  Yue Yu  Institute of Natural Sciences, School of Mathematical
Sciences, MOE-LSC, Shanghai Jiao Tong University, Shanghai, 200240, P. R.
China

###### Abstract

We present a simple new way–called Schrödingerisation– to simulate general
linear partial differential equations via quantum simulation. Using a simple
new transform, referred to as the warped phase transformation, any linear
partial differential equation can be recast into a system of Schrödinger’s
equations – in real time — in a straightforward way. This can be seen directly
on the level of the dynamical equations without more sophisticated methods.
This approach is not only applicable to PDEs for classical problems but also
those for quantum problems – like the preparation of quantum ground states,
Gibbs states and the simulation of quantum states in random media in the
semiclassical limit.

##  I Introduction

Quantum simulation is one of the most natural tasks for a quantum device — by
preparing outputs of Schrödinger’s equations directly via its own evolution
[1]. However, what can we say about the suitability of quantum systems in
simulating classical dynamical laws, typically in the form of ordinary or
partial differential equations (ODEs/PDEs) that are not Schrödinger’s
equations? This is of great importance for broader applications of quantum
computing in science and engineering.

While the most obvious applications are for classical problems, simulating
classical dynamics is also important for quantum problems. The most well-known
include the preparation of quantum ground states [2, 3] and Gibbs states [4],
which benefit from non-unitary dynamics, and are particularly important in
areas like quantum chemistry [5] and optimisation.

We know quantum simulation of classical dynamics is in principle possible, if
we accept reductionism and that the world is fundamentally quantum mechanical.
All classical dynamical laws are, with the exception of relativity, in
principle derivable from underlying Schrödinger’s equations. From a
computational perspective, we also know that any classical gate can be
trivially embedded into a quantum gate, so we expect quantum gates to be able
to simulate classical ones. However, the quantum degrees of freedom involved
could typically be significantly larger than that of the classical systems one
wishes to simulate. An essential question is then: are there less resource-
intensive ways to represent classical dynamics with Schrödinger’s equations?

Consider the case of numerical solutions to linear PDEs, which, upon spatial
and temporal discretizations, become a system of linear algebraic equations.
Quantum algorithm speedups in solving a system of linear algebraic equations
[6, 7] can subsequently lead to possible polynomial or super-polynomial
speedups in solving PDEs [8, 9]. For nonlinear counterparts see [10, 11].
However, here quantum simulation only plays the role of an algorithmic
primitive — due to discretisation in time, one does not prepare solution
states continuously in time $t$ by evolution of a Schrödinger equation in time
$t$.

To obviate the discretisation in time, a main difficulty is in finding a way
to represent non-unitary dynamics with a unitary one, which can be evolved by
quantum simulation like the Schrödinger equation. One way is via qubitisation
[12, 13] (or block-encoding) and this has very recently been applied to linear
PDEs [14]. This involves unitary dilation methods (e.g. [15]) and has origins
in quantum signal processing [16]. While this fairly general formalism can in
principle approximate the action of any non-unitary operator, it relies
heavily on building linear combinations of quantum states and finding suitable
polynomial approximations to the non-unitary operator, which is not always
simple to describe or to implement in practise.

For the heat equation, another method is the imaginary time evolution
approach, where the heat equation can be converted to a Schrödinger equation
through a change to imaginary time $t\rightarrow it$ [17]. However, the state
obeying Schrödinger’s equation and its imaginary time counterpart do not have
the same evolution. This means that, except for the steady state solution,
extra resources are necessary to map between the solution in the unitarily
evolving system to the other. For instance, one requires tomographic
measurements of quantum states at each small time step in the unitary
evolution [18]. Furthermore, the imaginary time evolution method does not
apply beyond the heat equation.

Thus we have an important question: Is there a simpler and generic way of
obtaining Schrödinger’s equations naturally from any given linear dynamics?

We propose a new paradigm — based on a simple transformation called the warped
phase transformation — that can map any linear PDE to Schrödinger equations in
real time. By the quantum simulation of this Schrödinger dynamics in time $t$,
it is possible to prepare solutions of the original PDEs at any time $t$, in
the form of a quantum state. We call this the Schrödingerisation approach.

To illustrate this approach, we present the heat equation example in detail.
With a warped phase transformation and a Fourier transform, the heat equation
can be transformed into a system of uncoupled Schrödinger equations with
rescaled time coordinates — in real time — thus providing a conceptual
alternative to imaginary time evolution methods. We show how the quantum
simulation of these Schrödinger’s equations can be used to prepare quantum
ground states and Gibbs states.

In another example, we apply the Schrödingerisation approach to a linear
radiative transport equation. This can be used to simulate quantum systems in
random media in the semi-classical limit.

Finally, we show how to generalise the Schrödingerisation approach to any
linear PDE.

##  II Background

Consider linear $(d+1)$-dimensional PDEs for $u=u(t,x)$ that is first-order in
time $t\geq 0$

$\displaystyle\partial_{t}u=\mathcal{L}u,\quad u(0,x)=u_{0}(x).$ (1)

Here $x\in\mathbb{R}^{d}$ is the position in $d$ dimensions,
$\partial_{(\cdot)}$ is the partial derivative with respect to $(\cdot)$ and
$\mathcal{L}$ is a linear differential operator with respect to $x$.

Schrödinger’s equation is such a linear PDE equation where $u$ is called the
wavefunction and $\mathcal{L}u=-i(-\nabla_{x}^{2}+V(x))u$, where
$\nabla_{x}^{2}$ is the Laplace operator with respect to $x$ and $V(x)$ is the
potential function. To find numerical solutions to this problem, one can
discretise in $x$, with uniform mesh sizes along each dimension $\Delta x=2/M$
where $M$ is a positive even integer. Let $\boldsymbol{u}(t)$ represent the
$M^{d}$-dimensional vector whose entries are $u(t,x_{i})$ at grid points with
$i=1,...,M^{d}$. By definition of the wavefunction, the $l_{2}$ norm
$\|\boldsymbol{u}(t)\|=1$. The entries to this vector can be considered as
amplitudes of the $M^{d}$-dimensional quantum state
$|u(t)\rangle=\sum_{i=1}^{M^{d}}u(t,x_{i})|i\rangle$ where $\\{|i\rangle\\}$
is an orthonormal basis set. When given a state spanned by
$\\{|0\rangle,|1\rangle\\}$ over the field $\mathbb{C}$, this is called a
qubit. Thus $|u(t)\rangle$ can be described by a system of $d\log_{2}(M)$
qubits. It is straightforward to see that $\boldsymbol{u}(t)$ satisfies the
linear system of ODEs

$\displaystyle i\frac{d\boldsymbol{u}(t)}{dt}=\boldsymbol{H}\boldsymbol{u}(t)$
(2)

where $\boldsymbol{H}$ is a $M^{d}\times M^{d}$ Hermitian matrix that results
from a discretisation of the Schrödinger Hamiltonian $-\nabla^{2}_{x}+V(x)$ in
$x$. The solution to Eq. (2) is clearly
$\boldsymbol{u}(t)=\exp(-i\boldsymbol{H}t)\boldsymbol{u}(0)$ and one has the
corresponding quantum state evolution
$|u(t)\rangle=\exp(-i\boldsymbol{H}t)|u(0)\rangle$. Quantum simulation
addresses the question of how to realise this $\exp(-i\boldsymbol{H}t)$
operator, when given $\boldsymbol{H}$, and an estimation of the resources
required. We focus on the time-independent setting for now.

Quantum simulation falls roughly into two categories: analogue and digital
[19]. In the analogue case, a quantum system that naturally realises
$\boldsymbol{H}$ or can be mapped to this Hamiltonian is found. However, if
such an analogue quantum system cannot be found, digital quantum simulation
methods still allows the simulation for more general $\boldsymbol{H}$. In this
case, one counts resources as the number of queries to some given black-boxes
and the number of two-qubit gates needed. This is referred to as the query and
gate complexities respectively.

Let $s$ be the sparsity of $\boldsymbol{H}$ (maximum number of non-zero
entries in each row) and $\|\boldsymbol{H}\|_{\text{max}}$ be its max-norm
(value of largest entry in absolute value). We denote the $(i,j)^{\text{th}}$
entry to $\boldsymbol{H}$ as $H_{ij}$. A common set of black-boxes used in
Hamiltonian simulation is known as the sparse access.

######  Definition 1.

Sparse access to Hermitian matrix $\boldsymbol{H}$ refers to two unitary
black-boxes $O_{M}$ and $O_{F}$ such that
$O_{M}|j\rangle|k\rangle|z\rangle=|j\rangle|k\rangle|z\oplus H_{jk}\rangle$
and $O_{F}|j\rangle|l\rangle=|j\rangle|F(j,l)\rangle$. Here the function $F$
takes the row index $j$ and a number $l=1,2,...,s$ and outputs the column
index of the $l^{\text{th}}$ non-zero elements in row $j$.

There are quantum simulation protocols in terms of query complexity that scale
linearly in $t$ [12] using sparse access or linearly in $t$ up to logarithmic
factors [20].

######  Lemma 2.

[20] Let $\tau=st\|\boldsymbol{H}\|_{\text{max}}$. Then
$\exp(-i\boldsymbol{H}t)$ acting on $m_{H}$ qubits can be simulated to within
error $\varepsilon$ has query complexity
$\mathcal{O}(\tau\log(\tau/\varepsilon)/(\log\log(\tau/\varepsilon)))$ with
gate complexity
$\mathcal{O}\Big{(}\tau(m_{H}+\log^{2.5}(\tau/\varepsilon))\log(\tau/\varepsilon)/(\log\log(\tau/\varepsilon))\Big{)}$.

Throughout the paper, we use $\tilde{\mathcal{O}}$ to denote $\mathcal{O}$
where logarithmic terms are ignored.

##  III Schrödingerisation of the heat equation

We begin with the initial value problem of the linear heat equation with a
source term

$\partial_{t}u=-Hu,\\\ \quad u(0,x)=u_{0}(x),$ (3)

where $H=-\nabla_{x}^{2}+V(x)$ and $t\geq 0$. We introduce a real one-
dimensional variable $p>0$ and define

$\displaystyle w(t,x,p)=e^{-p}u(t,x).$ (4)

This transformation – which we call the warped phase transformation – is a
crucial ingredient that will allow us to transform the heat equation to
Schrödinger equations by extending the heat equation solution to a higher
‘warped’ extra dimension. Such a transformation was used in [21] for a
completely different purpose: to develop efficient quantum computing
algorithms for uncertainty quantification problems in PDEs.

First observe that one can recover the solution to the heat equation using
$u(t,x)=\int_{0}^{\infty}w(t,x,p)\mathrm{d}p=\int_{-\infty}^{\infty}\chi(p)w(t,x,p)\mathrm{d}p$,
where $\chi(p)$ is the indicator function with $\chi(p)=1$ for $p\geq 0$ and
$\chi(p)=0$ for $p<0$. One can also recover $u$ via $u(t,x)=e^{p}W(t,x,p)$ for
any $p>0$. One can extend the domain of $p$ to $(-\infty,\infty)$, with
$w=w(t,x,p)$ satisfying $\partial_{t}w+(\nabla^{2}_{x}-V(x))\partial_{p}w=0$
with evenly extended initial condition $w(0,x,p)=\exp(-|p|u_{0}(x))$. Let
$\tilde{w}=\tilde{w}(t,x,\eta)$ be the Fourier transform of $w$ in $p$ and
$\eta\in\mathbb{R}$ be the Fourier mode. Then $\tilde{w}$ satisfies a system
of uncoupled Schrödinger equations

$\displaystyle i\partial_{t}\tilde{w}=\eta(-\nabla_{x}^{2}+V(x))\tilde{w},$
(5)

one for each $\eta$! See Appendix A for details. We call this the
Schrödingerised heat equation. The role of $\eta$ can be interpreted as a
recaling in time with $t\rightarrow t\eta$ for the Schrödinger equation
$i\partial_{t}\tilde{w}=(-\nabla^{2}_{x}+V(x))\tilde{w}$. The rescaled time
remains real-valued.

To solve these equations numerically, we discretise the system in $x$ and $p$,
but not in $t$. We choose uniform mesh sizes $\Delta x=2/M$ for the position
variable in each dimension, $\Delta p=2L/N$ for the $p$ variable,
$\Delta\eta=2L/N$ for the $\eta$ variable, where $M$ and $N$ are even positive
integers and $L>0$. We introduce the vector
$\boldsymbol{w}(t)=\boldsymbol{w}_{[x]}\otimes\boldsymbol{w}_{[p]}$ where the
elements of the vector $\boldsymbol{w}_{[x]}$ are $w(t,x_{i},p)$ with
$i=1,...,M^{d}$ labelling the $x$ grid points and the elements of
$\boldsymbol{w}_{[p]}$ are $w(t,x,p_{j})$ with $p_{j}=j\Delta p$,
$j=-N/2,...,N/2$. We define $\tilde{\boldsymbol{w}}=(\mathbf{1}^{\otimes
M^{d}}\otimes\mathcal{F}_{p})\boldsymbol{w}$ where $\mathcal{F}_{p}$ is the
discrete Fourier transform with respect to variable $p$. Then the
discretisation of the Schrödingerised heat equation in Eq. (5) becomes the
following system of $M^{d}N$ ODEs

$\displaystyle
i\frac{d}{dt}\tilde{\boldsymbol{w}}(t)=(\boldsymbol{H}\otimes\boldsymbol{D})\tilde{\boldsymbol{w}}(t)=\boldsymbol{H}_{\text{total}}\tilde{\boldsymbol{w}}(t).$
(6)

where $\boldsymbol{H}_{\text{total}}$ is a Hermitian matrix. Here
$\boldsymbol{H}=(-\boldsymbol{P}^{2}_{1}-...-\boldsymbol{P}^{2}_{d}+\boldsymbol{V})$
is also Hermitian, originating from the discretisation of the Schrödinger
Hamiltonian $H$ in Eq. (3). We define the matrices
$\boldsymbol{P}_{l}=\mathbf{1}^{\otimes l-1}\otimes
P_{l}\otimes\mathbf{1}^{\otimes d-l}$, $l=1,...,d$, where $P_{l}$ is the
discretisation of the momentum operator $-i\partial_{x}$ with respect to the
$l^{\text{th}}$ spatial variable. The matrix
$\boldsymbol{V}=\text{diag}(V(x_{1}),...,V(x_{M^{d}}))$ is a diagonal matrix
with entries $V(x_{i})$. The matrix
$\boldsymbol{D}=\text{diag}(\mu_{1},...,\mu_{N})$ is also a diagonal matrix,
with entries $\mu_{j}=\pi(j-N/2)$.

It is clear that
$\tilde{\boldsymbol{w}}(t)=\exp(-i\boldsymbol{H}_{\text{total}}t)\tilde{\boldsymbol{w}}(0)$
and $\boldsymbol{H}_{\text{total}}$, being Hermitian, can be interpreted as a
time-independent Hamiltonian. We can define the quantum state
$|\tilde{w}(t)\rangle=(1/\|\tilde{\boldsymbol{w}}(t)\|)\sum_{i=1}^{M^{d}}\sum_{j=-N/2}^{N/2}\tilde{w}(t,x_{i},\eta_{j})|i,j\rangle$
where
$\|\tilde{\boldsymbol{w}}(t)\|^{2}=\sum_{i=1}^{M^{d}}\sum_{j=-N/2}^{N/2}|\tilde{w}(t,x_{i},\eta_{j})|^{2}$
is the $l_{2}$ norm squared of $\tilde{\boldsymbol{w}}(t)$. From Eq. (6) one
can recover the quantum state with input quantum state $|\tilde{w}(0)\rangle$
through
$|\tilde{w}(t)\rangle=\exp(-i\boldsymbol{H}_{\text{total}}t)|\tilde{w}(0)\rangle$.
By applying an inverse quantum Fourier transform $\mathcal{F}_{p}^{-1}$, with
respect to $p$, onto the second register we obtain
$|w(t)\rangle=(\mathbf{1}^{\otimes
M^{d}}\otimes\mathcal{F}^{-1}_{p})|\tilde{w}(t)\rangle$, where
$|w(t)\rangle=(1/\|\boldsymbol{w}(t)\|)\sum_{i=1}^{M^{d}}\sum_{j=1}^{N}w(t,x_{i},p_{j})|i,j\rangle$
and $\|\boldsymbol{w}(t)\|=\|\tilde{\boldsymbol{w}}(t)\|$. From $|w(t)\rangle$
one can recover the quantum state of $u$ whose entries are proportional to the
solutions of the heat equation
$|u(t)\rangle=(1/\|\boldsymbol{u}(t)\|)\sum_{i=1}^{M^{d}}u(t,x_{i})|i\rangle$.
We can do this by either projecting $|w(t)\rangle$ onto
$\mathbf{1}\otimes\sum_{k=N/2}^{N}|k\rangle\langle k|$ (projecting only onto
$p>0$), or using amplitude amplification to boost the chance of retrieving
$|u(t)\rangle$ to probability
$\sim\|\boldsymbol{u}(0)\|/\|\boldsymbol{u}(t)\|$.

We call this the Schrödingerisation approach to prepare
$|u(t)\rangle\propto\exp(-\boldsymbol{H}t)|u(0)\rangle$, where $|u(t)\rangle$
evolves with respect to $\boldsymbol{H}$ instead of
$\boldsymbol{H}_{\text{total}}$. The above method is in fact applicable to any
Hermitian $\boldsymbol{H}$, with sparsity $s$ and max-norm
$\|\boldsymbol{H}\|_{\text{max}}$. We extend to non-Hermitian $\boldsymbol{H}$
in the last section.

######  Theorem 3.

Given sparse-access to a $D\times D$ Hermitian matrix $\boldsymbol{H}$ and the
unitary $U_{initial}$ that prepares the initial quantum state
$|u(0)\rangle=U_{initial}|0\rangle$. With the Schrödingerisation approach, the
state $|u(t)\rangle$ can be prepared to precision $\epsilon$ with query and
gate complexity
$\tilde{\mathcal{O}}((\|\boldsymbol{u}(0)\|/\|\boldsymbol{u}(t)\|)st\|\boldsymbol{H}\|_{\text{max}}/\epsilon)$.

###### Proof.

See Appendix B. ∎

###  III.1 Preparing quantum ground states and Gibbs states

Suppose one wants to prepare a $D$-dimensional ground state $|E_{0}\rangle$
from a given state $|u(0)\rangle=\sum_{j=0}^{D-1}\alpha_{j}|E_{j}\rangle$,
where $\alpha_{j}\in\mathbb{C}$, $\|\boldsymbol{u}(0)\|=1$ and
$\\{|E_{j}\rangle\\}$ are the non-degenerate orthonormal eigenstates of a
(Hermitian) $D\times D$ Hamiltonian $\boldsymbol{H}$. $\boldsymbol{u}(t)$ is
the vector whose entries are the amplitudes of $|u(t)\rangle$. If one evolves
$\boldsymbol{u}(t)=\exp(-\boldsymbol{H}t)\boldsymbol{u}(0)$ according to
$\partial_{t}\boldsymbol{u}=-\boldsymbol{H}\boldsymbol{u}$ with initial
condition $\boldsymbol{u}(0)$, then one can write
$|u(t)\rangle=(1/\|\boldsymbol{u}(t)\|)\sum_{j=0}^{D-1}\alpha_{j}\exp(-E_{j}t)|E_{j}\rangle$.
Then $|u(t)\rangle=\exp(-\boldsymbol{H}t)|u(0)\rangle$ where $\boldsymbol{H}$
is a $D\times D$ Hermitian matrix with sparsity $s$ and max-norm
$\|\boldsymbol{H}\|_{\text{max}}$. Assuming the presence of a non-zero
spectral gap $\Delta=E_{1}-E_{0}>0$, the convergence to the ground state
$|u(t)\rangle\rightarrow|E_{0}\rangle$ is exponentially fast. If
$t_{\text{final}}$ denotes the time where the quantum fidelity between the
target ground state $|E_{0}\rangle$ and the prepared state $|u(t)\rangle$ be
$F(|E_{0}\rangle,|u(t)\rangle)=1-\epsilon$ for a small $\epsilon>0$, then the
dominating contribution to $t_{\text{final}}$ is
$\mathcal{O}((1/\Delta)\ln(1/(\epsilon|\alpha_{0}|^{2})))$. Combining this
with Theorem 3, one sees the total query and gate complexity costs in
preparing the ground state to fidelity $1-\epsilon$ for $\epsilon\ll 1$ is
$\tilde{\mathcal{O}}(s\|\boldsymbol{H}\|_{\text{max}}/(|\alpha_{0}|\Delta\epsilon))$.
Here the scaling in $|\alpha_{0}|$ and $\Delta$ is comparable to non-heuristic
schemes like quantum phase estimation [2] and near-optimal lower bounds for
ground state preparation [3], up to logarithmic factors. The drawback is that
this scheme has an extra factor $1/\epsilon$ compared to the near-optimal
schemes. This originates from
$\|\boldsymbol{H}\otimes\boldsymbol{D}\|_{\text{max}}\sim\|\boldsymbol{H}\|_{\text{max}}/\epsilon$.
It will be interesting to see if improvements to this factor can be achieved.
For details see Appendix C.

To create the Gibbs state at temperature $T$ corresponding to the same
$D\times D$ Hamiltonian $\boldsymbol{H}$, we can use the Schrödingerisation
approach to prepare the normalised pure state
$|\Psi(\beta)\rangle=\sum_{k}\sqrt{\exp(-\beta E_{k})/Z}|E_{k}E_{k}\rangle$,
where $Z=\text{Tr}(\exp(-\beta\boldsymbol{H}))$ is the partition function for
$\boldsymbol{H}$ with $\beta=1/(k_{B}T)$ and $T$ is the temperature. Then by
tracing out one register we obtain the Gibbs state
$\rho_{Gibbs}(\beta)=\sum_{j=0}^{D-1}\exp(-\beta E_{j})|E_{j}\rangle\langle
E_{j}|=\text{Tr}_{1}(|\Psi(\beta)\rangle\langle\Psi(\beta)|)$. This quantum
simulation method differs from previous methods [4] that need to make use of
quantum phase estimation to prepare $|\Psi(\beta)\rangle$. We rewrite
$|\Psi(\beta)\rangle\propto\exp(-(\beta/2)(\boldsymbol{H}\otimes\mathbf{1})\sum_{j=0}^{D-1}|E_{j}E_{j}\rangle$
where we can identify $|u(t)\rangle=|\Psi(\beta)\rangle$ by the application of
$\exp(-(\boldsymbol{H}\otimes\mathbf{1})t)$, with time $t=\beta/2$, onto the
state $|u(0)\rangle=(1/\sqrt{D})\sum_{j=0}^{D-1}|E_{j}E_{j}\rangle$. The
latter we assume is given. This is equivalent to performing the simple
extension $\boldsymbol{H}\rightarrow\boldsymbol{H}\otimes\mathbf{1}$ in Eq.
(6), for any Hermitian $\boldsymbol{H}$. Now $\|u(t)\|=\sqrt{Z}$ and
$\|u(0)\|=1/\sqrt{D}$ where $D$ is the total number of energy levels. Then a
straightforward application of Theorem 3 means we can prepare
$\rho_{Gibbs}(\beta)$ to precision $\epsilon$ with query and gate complexity
$\tilde{\mathcal{O}}(s\|\boldsymbol{H}\|_{\text{max}}\beta\sqrt{D/Z}/\epsilon)$.
Heuristic methods aside, this coincides with the best-known scaling, to our
knowledge, with respect to $D$ and $Z$ [4, 13].

##  IV Schrödingerisation of linear transport equation

An interesting application of solving the linear transport equation (non-
unitary dynamics) for a quantum problem is in simulating quantum states in
random media in the semi-classical regime. Here we are interested in the
evolution of a quantum system – described by the Schrödinger equation – in the
presence of a random potential $U(x)=U_{0}(x)+U_{1}(x/\zeta)$. Here $U_{0}(x)$
is a deterministic, slowly moving background potential and $U_{1}$ is the
random fluctuation depending on $x/\zeta$, the fast spatial variable, and
$R(|x-y|)=\langle U_{1}(x)U_{2}(y)\rangle$ is the covariance of the
fluctuations. We assume homogeneity and isotropy of the medium, where the
differential scattering cross-section $\sigma(k,k^{\prime})\propto
R(k-k^{\prime})\delta(k^{2}-(k^{\prime})^{2})$ with $\tilde{R}(|k|)$ being the
power spectrum and $\Sigma(k)$ denotes the total scattering cross-section. In
the semiclassical limit, the Wigner transform of the wave function will
converge to the particle probability density function $W=W(t,x,k)$ governed by
the linear transport equation [22]

$\displaystyle\partial_{t}W+k\cdot\nabla_{x}W=\int
dk^{\prime}\sigma(k,k^{\prime})W(t,x,k^{\prime})-\Sigma(k)W.$ (7)

By Schrödingerising this equation, we can prepare, to precision $\epsilon$,
the quantum state $|W(t)\rangle$ whose amplitudes are proportional
$\boldsymbol{W}(t)$, which is the vector formed by taking $W$ at mesh points
$(x_{i},k_{j})$ formed by discretising $x$ and $k$. Following the
Schrödingerisation approach (see Appendix D for details), the corresponding
$\boldsymbol{H}_{\text{total}}$ for this problem is
$\boldsymbol{H}_{\text{total}}=\boldsymbol{L}_{[\xi,k]}\otimes\boldsymbol{1}-\mathbf{1}\otimes\boldsymbol{\Sigma}\otimes\boldsymbol{D}+\mathbf{1}\otimes\boldsymbol{\sigma}\otimes\boldsymbol{D}$,
where $\boldsymbol{\sigma}$ and $\boldsymbol{\Sigma}$ are the differential and
total cross-section matrices. Matrix $\boldsymbol{L}_{[\xi,k]}$ has entries
corresponding to $\xi\cdot k$, where $\xi$ is the Fourier mode of $x$ and
arises from discretising $k\cdot\nabla_{x}W$, which is a fundamental part of
the transfer equation in Eq. (7). This part exists whether or not the
Schrödingerisation approach is used. Here
$\|\boldsymbol{L}_{[\xi,k]}\|_{\text{max}}\sim\|\boldsymbol{\sigma}\otimes\boldsymbol{D}\|_{\text{max}}\sim\mathcal{O}(1/\epsilon)$
where $\boldsymbol{\sigma}$ and $\boldsymbol{\Sigma}$ are order one factors.
Since the warped transformation only adds one derivative in $p$ in the
scattering term, and the new transformed equation remains first order, then
the order of $\|\boldsymbol{H}_{\text{total}}\|_{\text{max}}$ remains the same
order in $1/\epsilon$ as the original transport equation. This means that the
Schrödingerisation approach does not give rise to an extra $O(1/\epsilon)$
factor compared to the original transport equation. Then to retrieve
$|W(t)\rangle$, it can be shown that the query and gate complexities are
respectively
$\sim\tilde{\mathcal{O}}((\|\boldsymbol{W}(0)\|/\|\boldsymbol{W}(t)\|)s(\boldsymbol{\sigma})/\epsilon)$
and
$\sim\tilde{\mathcal{O}}((\|\boldsymbol{W}(0)\|/\|\boldsymbol{W}(t)\|)ds(\boldsymbol{\sigma})/\epsilon)$,
where $s(\boldsymbol{\sigma})$ is the sparsity of the differential cross-
section matrix. This approach can also be applied to find the stationary state
for $|W(t)\rangle$.

Here $d\geq 3$ and refers to the spatial dimension for Schrödinger’s equation.
It is also possible to interpret this in terms of multiple $n=d/3$ particles
each moving in three spatial dimensions so a high $d$ limit is a large $n$
limit. See Appendix D.

While $|W(t)\rangle$ not the usual output of a quantum simulator, which
outputs directly a wavefunction $|\psi(t)\rangle$, it is still possible to
recover the semiclassical limit of the physical observables of
$|\psi(t)\rangle$ by taking the moments of $W$. Namely, any expectation value
$\langle\hat{G}\rangle=\langle\psi(t)|\hat{G}|\psi(t)\rangle=\iint
dxdkW(t,x,k)g(x,k)$, where $g(x,k)=1,k,k^{2}/2$, gives mass, momentum and
energy respectively. Given access to the quantum state $|g(x,k)\rangle$ whose
amplitudes are proportional to $g(x,k)$ at mesh points $(x_{i},k_{j})$,
$\langle\hat{G}\rangle$ is then proportional to the quantum fidelity between
$|g(x,k)\rangle$ and $|W(t,x,k)\rangle$. This fidelity can be recovered in
different ways, for instance, via a quantum swap-test [23].

##  V General formulation

The Schrödingerisation approach can in fact be generalised to any linear
$(d+1)$-dimensional PDE for $u(t,x)$. We study the PDE that is first-order in
time. A system with higher-order derivatives in $t$ can be written as an
enlarged system with first-order time-derivatives by introducing a new
variable $u^{\prime}=\partial_{t}u$. Using the same discretisation scheme as
before, where the entries of vector $\boldsymbol{u}(t)$ with size $M^{d}$ are
solutions to a homogeneous linear PDE on the grid points $(t,x_{i})$, we are
left with a system of linear ODEs

$\displaystyle d\boldsymbol{u}(t)/dt=-\boldsymbol{A}\boldsymbol{u}(t),$ (8)

where the matrix $\boldsymbol{A}$ is in general not Hermitian, unlike in the
heat equation example. However, it is always possible to make the
decomposition $\boldsymbol{A}=\boldsymbol{H}+i\boldsymbol{\bar{H}}$ in terms
of Hermitian matrices
$\boldsymbol{H}=(\boldsymbol{A}+\boldsymbol{A}^{\dagger})/2$ and
$\boldsymbol{\bar{H}}=i(\boldsymbol{A}^{\dagger}-\boldsymbol{A})/2$. To ensure
stability, we can assume $\boldsymbol{H}$ to be positive semi-definite. The
heat equation example belongs to the case $\boldsymbol{\bar{H}}=0$, whereas
purely Schrödinger dynamics corresponds to $\boldsymbol{H}=0$. A
generalisation to inhomogeneous linear PDEs is possible by a straightforward
dilation of $\boldsymbol{u}$ and $\boldsymbol{A}$. This is also useful for
solving boundary value problems [24] by quantum simulation.

We can apply our warped phase transformation here to $\boldsymbol{u}(t)$ in
the $p>0$ region, $\boldsymbol{v}(t,p)=e^{-p}\boldsymbol{u}(t)$, where the
original solution can be recovered by using
$\boldsymbol{u}(t)=\int_{-\infty}^{\infty}\chi(p)\boldsymbol{v}(t,p)dp$. Just
like for the heat equation, it is possible to extend the initial condition to
$p<0$ by defining $\boldsymbol{v}(0,p)=\exp(-|p|)\boldsymbol{u}(0)$ and in the
$p\in(-\infty,\infty)$ region $\boldsymbol{v}(t,p)$ satisfies
$\partial_{t}\boldsymbol{v}+\boldsymbol{H}\partial_{p}\boldsymbol{v}-i\boldsymbol{\bar{H}}\boldsymbol{v}=0$.
Let $\tilde{\boldsymbol{v}}(t,\eta)$, $\eta\in\mathbb{R}$ be the Fourier
transform of $\boldsymbol{v}$ in $p$, we obtain the generalisation of Eq. (6)

$\displaystyle
i\partial_{t}\tilde{\boldsymbol{v}}=(\eta\boldsymbol{H}+\boldsymbol{\bar{H}})\tilde{\boldsymbol{v}}$
(9)

which is still a system of Schrödinger equations, one for each $\eta$, since
$\eta\boldsymbol{H}+\boldsymbol{\bar{H}}$ is Hermitian. Note this form is not,
in a strict sense, Schrödinger’s equations per se, where the Hamiltonian is
usually a Laplacian plus a potential function. However, the first order
differential operator can be cast explicitly into the Schrödinger form by
making another warped phase transformation. We omit the details here since
it’s not crucial for quantum simulation to be possible. As previously, we can
proceed by discretising $\eta$ to obtain

$\displaystyle
i\frac{d}{dt}\tilde{\boldsymbol{v}}=(\boldsymbol{H}\otimes\boldsymbol{D}+\boldsymbol{\bar{H}}\otimes\mathbf{1})\tilde{\boldsymbol{v}}=\boldsymbol{H}_{\text{total}}\tilde{\boldsymbol{v}},$
(10)

and the state evolves
$\tilde{\boldsymbol{v}}(t)=\exp(i\boldsymbol{H}_{\text{total}}t)\tilde{\boldsymbol{v}}(0)$
with respect to unitary dynamics. Let
$s\sim\max(s(\boldsymbol{H}),s({\boldsymbol{\bar{H}}}))$ be the maximum of the
sparsity of $\boldsymbol{H}$ and $\boldsymbol{\bar{H}}$. The max-norm
$\|\boldsymbol{H}_{\text{total}}\|_{\text{max}}\sim\max(\|\boldsymbol{H}\|_{\text{max}}/\epsilon,\|\boldsymbol{\bar{H}}\|_{\text{max}})$.
For instance, in the transport equation case, the two terms in the latter
expression are of the same order. Then we can simulate $|u(t)\rangle$ with the
following resources, where the proof follows in the same way as Theorem 3 with
$D=M^{d}$.

######  Theorem 4.

Given sparse access to the $M^{d}\times M^{d}$ matrix
$\boldsymbol{H}_{\text{total}}$ and the unitary $U_{initial}$ that prepares
the initial quantum state $|\boldsymbol{u}(0)\rangle$ to precision $\epsilon$.
With the Schrödingerisation approach, the state $|\boldsymbol{u}(t)\rangle$
can be prepared with query complexity
$\tilde{\mathcal{O}}((\|\boldsymbol{u}(0)\|/\|\boldsymbol{u}(t)\|)st\|\boldsymbol{H}_{\text{total}}\|_{\text{max}})$
and
$\tilde{\mathcal{O}}(\|\boldsymbol{u}(0)\|/\|\boldsymbol{u}(t)\|dts\|\boldsymbol{H}_{\text{total}}\|_{\text{max}})$
additional two-qubit gates.

##  VI Discussion

We have introduced a conceptually new method, called Schrödingerisation, which
makes it possible to simulate solutions of any linear partial differential
equations using quantum simulation. This method is framed in the traditional
language of dynamical equations. The connection between classical dynamics and
its corresponding Schrödinger’s equations can be seen more easily at the level
of the dynamical equations without more sophisticated methods.

For more in-depth technical details, discussions and applications to other
linear PDEs like Liouville, Fokker-Planck, Vlaslov-Fokker-Planck, Black-
Scholes equations and nonlinear ODEs, please see our technical companion paper
[24].

##  VII Acknowledgements

SJ was partially supported by the NSFC grant No. 12031013, the Shanghai
Municipal Science and Technology Major Project (2021SHZDZX0102), and the
Innovation Program of Shanghai Municipal Education Commission (No.
2021-01-07-00-02-E00087). NL acknowledges funding from the Science and
Technology Program of Shanghai, China (21JC1402900).

## References

  * Feynman [2018] R. P. Feynman, Simulating physics with computers, in  _Feynman and computation_ (CRC Press, 2018) pp. 133–153.
  * Ge  _et al._ [2019] Y. Ge, J. Tura, and J. I. Cirac, Faster ground state preparation and high-precision ground energy estimation with fewer qubits, Journal of Mathematical Physics 60, 022202 (2019).
  * Lin and Tong [2020] L. Lin and Y. Tong, Near-optimal ground state preparation, Quantum 4, 372 (2020).
  * Poulin and Wocjan [2009] D. Poulin and P. Wocjan, Sampling from the thermal quantum gibbs state and evaluating partition functions with a quantum computer, Physical review letters 103, 220502 (2009).
  * Lloyd [1996] S. Lloyd, Universal quantum simulators, Science 273, 1073 (1996).
  * Harrow  _et al._ [2009] A. W. Harrow, A. Hassidim, and S. Lloyd, Quantum algorithm for linear systems of equations, Physical review letters 103, 150502 (2009).
  * Childs  _et al._ [2017] A. M. Childs, R. Kothari, and R. D. Somma, Quantum algorithm for systems of linear equations with exponentially improved dependence on precision, SIAM Journal on Computing 46, 1920 (2017).
  * Childs  _et al._ [2021] A. M. Childs, J.-P. Liu, and A. Ostrander, High-precision quantum algorithms for partial differential equations, Quantum 5, 574 (2021).
  * Montanaro and Pallister [2016] A. Montanaro and S. Pallister, Quantum algorithms and the finite element method, Physical Review A 93, 032324 (2016).
  * Liu  _et al._ [2021] J.-P. Liu, H. Ø. Kolden, H. K. Krovi, N. F. Loureiro, K. Trivisa, and A. M. Childs, Efficient quantum algorithm for dissipative nonlinear differential equations, Proceedings of the National Academy of Sciences 118, e2026805118 (2021).
  * Jin and Liu [2022] S. Jin and N. Liu, Quantum algorithms for computing observables of nonlinear partial differential equations, arXiv:2202.07834 (2022).
  * Low and Chuang [2019] G. H. Low and I. L. Chuang, Hamiltonian simulation by qubitization, Quantum 3, 163 (2019).
  * Gilyén  _et al._ [2019] A. Gilyén, Y. Su, G. H. Low, and N. Wiebe, Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics, in  _Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing_ (2019) pp. 193–204.
  * An  _et al._ [2022] D. An, J.-P. Liu, D. Wang, and Q. Zhao, A theory of quantum differential equation solvers: limitations and fast-forwarding, arXiv:2211.05246 (2022).
  * Horn and Johnson [2012] R. A. Horn and C. R. Johnson,  _Matrix analysis_ (Cambridge university press, 2012).
  * Low  _et al._ [2016] G. H. Low, T. J. Yoder, and I. L. Chuang, Methodology of resonant equiangular composite quantum gates, Physical Review X 6, 041067 (2016).
  * Lehtovaara  _et al._ [2007] L. Lehtovaara, J. Toivanen, and J. Eloranta, Solution of time-independent schrödinger equation by the imaginary time propagation method, Journal of Computational Physics 221, 148 (2007).
  * Motta  _et al._ [2020] M. Motta, C. Sun, A. T. Tan, M. J. O’Rourke, E. Ye, A. J. Minnich, F. G. Brandão, and G. K. Chan, Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution, Nature Physics 16, 205 (2020).
  * Daley  _et al._ [2022] A. J. Daley, I. Bloch, C. Kokail, S. Flannigan, N. Pearson, M. Troyer, and P. Zoller, Practical quantum advantage in quantum simulation, Nature 607, 667 (2022).
  * Berry  _et al._ [2015] D. W. Berry, A. M. Childs, and R. Kothari, Hamiltonian simulation with nearly optimal dependence on all parameters, in  _2015 IEEE 56th annual symposium on foundations of computer science_ (IEEE, 2015) pp. 792–809.
  * Golse  _et al._ [2022] F. Golse, S. Jin, and N. Liu, Quantum algorithms for uncertainty quantification: application to partial differential equations, arXiv:2209.11220 (2022).
  * Ryzhik  _et al._ [1996] L. Ryzhik, G. Papanicolaou, and J. B. Keller, Transport equations for elastic and other waves in random media, Wave Motion 24, 327 (1996).
  * Buhrman  _et al._ [2001] H. Buhrman, R. Cleve, J. Watrous, and R. De Wolf, Quantum fingerprinting, Physical Review Letters 87, 167902 (2001).
  * Jin  _et al._ [2022] S. Jin, N. Liu, and Y. Yu, Quantum simulation of partial differential equations via schrödingerisation: technical details, upcoming arXiv XXXX (2022).
  * Brassard  _et al._ [2002] G. Brassard, P. Hoyer, M. Mosca, and A. Tapp, Quantum amplitude amplification and estimation, Contemporary Mathematics 305, 53 (2002).
  * Erdős and Yau [2000] L. Erdős and H.-T. Yau, Linear boltzmann equation as the weak coupling limit of a random schrödinger equation, Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences 53, 667 (2000).
  * Mouhot and Neumann [2006] C. Mouhot and L. Neumann, Quantitative perturbative study of convergence to equilibrium for collisional kinetic models in the torus, Nonlinearity 19, 969 (2006).
  * Chandrasekhar [1960] S. Chandrasekhar,  _Radiative transfer_ (Dover Publications, Inc., New York, 1960) pp. xiv+393.

##  Appendix A Appendix A: Warped phase transformation

Here we define in the $p>0$ region

$w(t,x,p)=\mathrm{e}^{-p}u(t,x),\qquad p>0.$

A simple calculation shows

$\partial_{t}w+\partial_{p}\nabla_{x}w=0,\qquad p>0.\\\ $ (11)

We can recover $u$ using $w$ via

$u(t,x)=\int_{0}^{\infty}w(t,x,p)\mathrm{d}p=\int_{-\infty}^{\infty}\chi(p)w(t,x,p)\mathrm{d}p,$
(12)

where $\chi(p)=1$ for $p>0$ and $\chi(p)=0$ for $p<0$. Alternatively, one can
also choose any $p_{*}>0$ and let

$u(t,x)=\mathrm{e}^{p_{*}}w(t,x,p_{*}).$ (13)

Now we want to see that no boundary condition is required at $p=0$. Assume $x$
is defined in a periodic domain and we apply the Fourier transform on $w$ with
respect to $x$, denoted $\hat{w}(t,\xi,p)$, with Fourier modes
$\xi=(\xi_{1},...,\xi_{d})^{T}$. One then arrives at a convection equation
$\partial_{t}\hat{w}-|\xi|^{2}\partial_{p}\hat{w}=0$, where
$|\xi|^{2}=\xi_{1}^{2}+\cdots+\xi_{d}^{2}$. Here the solution $\hat{w}$ moves
from the right to the left, so no boundary condition is needed at $p=0$.

For the convenience of numerical approximation, we extend the domain to $p<0$.
This extension will not affect the solution of $w$ in the region $p>0$ since
$w$ convects from right to the left. Thus, we can symmetrically extend the
initial data of $w$ to $p<0$ and keep Eq. (11):

$\begin{cases}\partial_{t}w+\Delta_{x}\partial_{p}w=0,\qquad
p\in(-\infty,\infty),\\\ w(0,x,p)=\mathrm{e}^{-|p|}u_{0}(x).\end{cases}$ (14)

The only difference is that now the initial condition is
$w(0,x,p)=\exp(-|p|)u_{0}(x)$ for $p\in(-\infty,\infty)$ instead of just
$w(0,x,p)=\exp(-p)u_{0}(x)$ for $p>0$.

The solution coincides with the solution of Eq. (11) when $p>0$. Due to the
exponential decay of $\mathrm{e}^{-|p|}$ one can (computationally) impose the
periodic boundary condition $w(t,x,p=-L)=w(t,x,p=L)=0$ along the $p$-direction
for some $L>0$ suitably large but finite. Then the Fourier transform on $w$
with respect to $p$ gives

$\partial_{t}\tilde{w}-\mathrm{i}\eta\nabla^{2}_{x}\tilde{w}=0\qquad\mbox{or}\qquad\mathrm{i}\partial_{t}\tilde{w}=-\eta\nabla^{2}_{x}\tilde{w},$
(15)

where $\tilde{w}(t,x,\eta)$, $\eta\in\mathbb{R}$, is the Fourier transform of
$w$ in $p$. Equation (15) is clearly the Schrödinger equation, for every
$\eta$!

For more details and numerical experiments on justifying the setup in Eq.
(14), see our technical companion paper [24].

##  Appendix B Appendix B: Proof of Theorem 3

From Eq. (6) one can prepare the $\tilde{\mathcal{O}}(\log(D))$-qubit state
$|\tilde{w}(t)\rangle=\exp(-i\boldsymbol{H}_{\text{total}}t)|\tilde{w}(0)\rangle$
by simulating the Hamiltonian
$\boldsymbol{H}_{\text{total}}=\boldsymbol{H}\otimes\boldsymbol{D}$ to time
$t$. Here the results hold so long as $\boldsymbol{H}$ is Hermitian. Let $s$
and $\|\boldsymbol{H}\|_{\text{max}}$ be the sparsity and max-norm of a more
general $D\times D$ Hermitian matrix $\boldsymbol{H}$.

Then, given $|\tilde{w}(0)\rangle$, one can apply Lemma 2, so that one
prepares $|\tilde{w}(t)\rangle$ with query and gate complexity
$\tilde{\mathcal{O}}(s\|\boldsymbol{H}_{\text{total}}\|_{\text{max}}t)$, where
$\tilde{\mathcal{O}}$ denotes that logarithmic terms, like $\log(D)$ and
$\log(1/\epsilon)$, are ignored. Now
$\|\boldsymbol{H}_{\text{total}}\|_{\text{max}}=\|\boldsymbol{H}\|_{\text{max}}\|\boldsymbol{D}\|_{\text{max}}$,
where $\|\boldsymbol{D}\|_{\text{max}}=\mathcal{O}(N)$ by definition and
$N=\mathcal{O}(1/\epsilon)$. Thus
$\tilde{\mathcal{O}}(st\|\boldsymbol{H}_{\text{total}}\|_{\text{max}})=\tilde{\mathcal{O}}(st\|\boldsymbol{H}\|_{\text{max}}/\epsilon)$.

From $|\tilde{w}(t)\rangle$ that results from the Hamiltonian simulation
problem, the simplest way to prepare $|u(t)\rangle$ from
$|\tilde{w}(t)\rangle$ is to project $|\tilde{w}(t)\rangle$ with operator
$\hat{P}=\mathbf{1}\otimes\sum_{k=N/2}^{N}|k\rangle\langle k|$ (which projects
only onto $p>0$ states). This leads to

$\displaystyle\hat{P}|w(t)\rangle=\frac{\|u(t)\|\|\exp(-p)\|}{\|w(t)\|}|\exp(-p)\rangle|u(t)\rangle$
(16)

where $|\exp(-p)\rangle=(1/\|\exp(-p)\|)\sum_{k=N/2}^{N}\exp(-p_{k})|k\rangle$
with normalisation constant
$\|\exp(-p)\|^{2}=\sum_{k=N/2}^{N}\exp(-2p_{k})\sim\mathcal{O}(N\int_{0}^{\infty}\exp(-2p)dp)=\mathcal{O}(N)$.
Thus a simple projection retrieves $|u(t)\rangle$ with probability
$(\|u(t)\|\|\exp(-p)\|/\|w(t)\|)^{2}\sim N(\|u(t)\|/\|w(t)\|)^{2}$. We also
used $\|\tilde{w}(t)\|=\|w(t)\|$.

One can also use amplitude amplification to boost this probability to
$\sqrt{N}\|u(t)\|/\|w(t)\|$ in the usual manner, by assuming access to the
oracle $Q=-S_{w}S_{p}$ where $S_{w}=\mathbf{1}-|w(t)\rangle\langle w(t)|$ and
$S_{p}=\mathbf{1}-2\hat{P}$ [25]. This means
$\tilde{\mathcal{O}}(\|w(t)\|/(\sqrt{N}\|u(t)\|))$ queries to $Q$ is
sufficient. Combining this with Lemma 2, we see that this requires a total
query complexity of
$\tilde{\mathcal{O}}(st\|\boldsymbol{H}\|_{\text{max}}\|w(t)\|/(\sqrt{N}\|u(t)\|))$.

It can be shown, using the quadrature rule, that
$(1/(M^{d}N))\|w(t)\|^{2}\approx\iint_{-\infty}^{\infty}dxdp|w(t,x,p)|^{2}=\int_{-\infty}^{\infty}|u(0)|^{2}\approx(1/M^{d})\|u(0)\|^{2}$.
Thus $\|w(t)\|/(\sqrt{N}\|u(t)\|)=\|u(0)\|/\|u(t)\|$.

In the above, we assumed access to $|\tilde{w}(0)\rangle$, while it is more
natural to assume access to $|u(0)\rangle$. It is straightforward to obtain
the initial state $|\tilde{w}(0)\rangle$ from $|u(0)\rangle$. Since
$w(0,x,p)=\exp(-|p|)u(0,x)$, then $|w(0)\rangle=|u(0)\rangle|\exp(-p)\rangle$,
where we assume we can prepare the state $|\exp(-p)\rangle$. Then applying the
quantum Fourier transform to $|w(0)\rangle$, with respect to $p$, gives
$|\tilde{w}(0)\rangle$. To perform an inverse Fourier transform, only
$\tilde{\mathcal{O}}(1)$ two-qubit gates are required, since we have only
$\log(N)$ qubits in the $p$-variable. The number of qubits in
$|\tilde{w}(t)\rangle$ is $\log(DN)$.

Thus, when given $|u(0)\rangle$, the total query and gate complexity to
prepare $|u(t)\rangle$ is
$\tilde{\mathcal{O}}((\|u(0)\|/\|u(t)\|)st\|\boldsymbol{H}\|_{\text{max}}/\epsilon)$.

##  Appendix C Appendix C: Preparation of quantum ground state via
Schrödingerisation

Given a (Hermitian) Hamiltonian $\boldsymbol{H}$, an important question is how
to prepare its ground state. If given a quantum state $|u(0)\rangle$, where
$\|u(0)\|=1$, then unitary evolution with respect to a time-independent
$\boldsymbol{H}$ for time $\tau$ gives
$|u(\tau)\rangle=\exp(-i\boldsymbol{H}\tau)|u(0)\rangle$. Here it is not clear
at what time $t$ one might be able to retrieve the ground state, if at all. So
an extension to non-unitary evolution would be helpful.

One method is the imaginary time evolution approach, which maps
$\tau\rightarrow t=i\tau$. In this case, one can write the non-unitary
evolution as

$\displaystyle|u(t)\rangle=\frac{e^{-\boldsymbol{H}t}}{\|u(t)\|}|u(0)\rangle.$
(17)

Since this is now non-unitary evolution $\|u(t)\|\neq\|u(t=0)\|$. However,
there is now no direct way to transform the solution $|u(\tau)\rangle$ into
Eq. (17) since there is no $\tau\rightarrow t$ conversion that can be used
once the solutions are obtained. The two solutions will only agree in the
stationary limit $\partial u/\partial t=0$. One can instead use an indirect
method [18], where for each $t$, one finds the corresponding Hamiltonian
evolving with $\tau$, by making tomographic measurements of $|u\rangle$ at
every small time-step. We seek a more direct method, without intervening the
quantum system during the process, that still allows us to perform this non-
unitary evolution.

What we lose in unitarity of evolution in Eq. (17), we gain in the simplicity
of seeing ground states emerge. Let the orthonormal eigenfunctions of
$\boldsymbol{H}$ be $\\{E_{j}\rangle\\}$ for $j=1,...,D$. We can thus
decompose $|u(0)\rangle$

$\displaystyle|u(0)\rangle=\sum_{j=0}^{D-1}\alpha_{j}|E_{j}\rangle,\quad\alpha_{j}\in\mathbb{C}$
(18)

with corresponding eigenvalues $E_{j}$ where
$\sum_{j=0}^{D-1}|\alpha_{j}|^{2}=1$. Then it is clear, if we evolve non-
unitarily,

$\displaystyle|u(t)\rangle=\frac{e^{-\boldsymbol{H}t}}{\|u(t)\|}|u(0)\rangle=\frac{1}{\|u(t)\|)}\sum_{j=0}^{D-1}\alpha_{j}e^{-E_{j}t}|E_{j}\rangle$
(19)

This means that for low eigenenergies, the decay is the slowest and
$|E_{0}\rangle$ would emerge as $t\rightarrow\infty$ for gapped Hamiltonians.
Furthermore, the convergence rate toward the ground state is exponential with
the decay rate of $\exp((E_{0}-E_{1})t)$, if there is a spectral gap, i.e.,
$E_{1}>E_{0}$.

Here we can use Schrödingerisation to obtain the state
$|u(t)\rangle\propto\exp(-\boldsymbol{H}t)|u(0)\rangle$. To specify how long
$t$ we need to evolve, let $t_{\text{final}}$ be the minimum evolution time so
that quantum fidelity between the prepared $|u(t_{\text{final}})\rangle$ state
and the true ground state $|E_{0}\rangle$ is greater or equal to $1-\epsilon$,
i.e.,

$\displaystyle|\langle
u(t_{\text{final}})|E_{0}\rangle|^{2}=\frac{|\alpha_{0}|^{2}e^{-2E_{0}t_{\text{final}}}}{\|u(t_{\text{final}})\|^{2}}\geq
1-\epsilon$ (20)

for some $\epsilon\ll 1$. For a normalised initial state $\|u(0)\|=1$ we have

$\displaystyle\|u(t)\|=\|e^{-\boldsymbol{H}t}|u(0)\rangle\|.$ (21)

thus
$\|u(t)\|^{2}=|\alpha_{0}|^{2}\exp(-2E_{0}t)(1+(|\alpha_{1}|/|\alpha_{0}|)^{2}\exp(-2\Delta
t))+L$ where $0<\Delta=E_{1}-E_{0}$ is the spectral gap and we can assume
$L=\sum_{k\geq 2}|\alpha_{k}|^{2}e^{-2tE_{k}}\ll 1$ for $t\gg 1$. Then for
$\delta\ll 1$, one can insert this into Eq. (20) and obtain

$\displaystyle
t_{\text{final}}\gtrsim(1/2\Delta)\ln(|\alpha_{1}|^{2}/(\epsilon|\alpha_{0}|^{2})\sim(1/\Delta)\ln(1/(\epsilon|\alpha_{0}|^{2}))$
(22)

when $|\alpha_{k}|\ll 1$ for all $k\geq 2$.

Inserting $t=t_{\text{final}}$ in Theorem 3 and $E_{0}\ll t_{\text{final}}$,
one gets the total query and gate complexity
$\tilde{\mathcal{O}}(s\|\boldsymbol{H}\|_{\text{max}}/(|\alpha_{0}|\Delta\epsilon))$
for preparing the ground state of $\boldsymbol{H}$.

The proof for Gibbs state preparation proceeds similarly, except there
$t\propto 1/T$ where $T$ is the temperature of the Gibbs state.

##  Appendix D Appendix D: Transport equation and quantum systems in random
media

Energy conservation for general wave problems is given by the transport
equation for the scalar energy density

$\displaystyle\partial_{t}f+\nabla_{k}\omega\cdot\nabla_{x}f-\nabla_{x}\omega\cdot\nabla_{k}f=S,$
(23)

with initial condition $f(0,x,k)=f_{0}(x,k)$, where $f=f(t,x,k)$ is the scalar
energy density and $\omega=\omega(x,k)$ is the frequency of the wave with wave
vector $k\in\mathbb{S}^{d-1}$ at point $x\in[-1,1]^{d}$. We assume periodic
boundary conditions are imposed. Here $S=\int
dk^{\prime}(\sigma(x,k,k^{\prime})f(t,x,k^{\prime})-\Sigma(x,k)f)$ where
$\sigma=\sigma(x,k,k^{\prime})$ is the differential scattering cross section
and $\Sigma(x,k)=\int dk^{\prime}\sigma(x,k^{\prime},k)$ is the total cross-
section.

An interesting example of this equation for a quantum problem appears in the
simulation of quantum states in random media in the semi-classical regime.
Here we study the evolution of a quantum system in the presence of a random
potential $U(x)=U_{0}(x)+U_{1}(x/\zeta)$ where $U_{0}(x)$ is a deterministic
slowly moving background potential, $U_{1}$ is the random fluctuation
depending on $x/\zeta$, the fast spatial variable.

Let $R(|x-y|)=\langle U_{1}(x)U_{2}(y)\rangle$ be the covariance of the
fluctuations, which are assumed to be spatially homogeneous and isotropic, and
$R(|k|)=(1/(2\pi))^{d}\int dk\exp(iky)R(x)$ is the power spectrum. The
differential cross-section and the total cross-section in this case are
respectively
$\sigma(k,k^{\prime})/(4\pi)=R(k-k^{\prime})\delta(k^{2}-(k^{\prime})^{2})$
and $\Sigma(k)/(4\pi)=\int
dk^{\prime}R(k-k^{\prime})\delta(k^{2}-(k^{\prime})^{2})$. For the Schrödinger
equation, $\omega^{2}=k^{2}/2+U_{0}(x)=k^{2}/2$, we can ignore the slowly-
moving background potential $U_{0}(x)$.

The semiclassical or weakly coupling limit of the Wigner function for the
quantum system in this random medium, which is the particle density
distribution, solves the transport equation in Eq. (23) [22] and for any
$d\geq 3$ the above forms of the cross-sections hold (Born approximation) [26]
and

$\displaystyle\partial_{t}W+k\cdot\nabla_{x}W=\int
dk^{\prime}\sigma(k,k^{\prime})W(t,x,k^{\prime})-\Sigma(k)W.$ (24)

We can apply the warped phase transformation $V=V(t,x,k,p)=\exp(-p)W$ to get

$\displaystyle\partial_{t}V+k\cdot\nabla_{x}V=-\int
dk^{\prime}\sigma(k,k^{\prime})\partial_{p}V(t,x,k^{\prime},p)+\Sigma(k)\partial_{p}V.$
(25)

Now apply a Fourier transform with respect to both the $x$ and $p$ variables,
which have respective conjugate variables $\xi\in\mathbb{R}^{d}$ and
$\eta\in\mathbb{R}$ and define this Fourier transform as

$\displaystyle\tilde{V}=\tilde{V}(t,\xi,k,\eta)=\int dxdp\,e^{-i\eta
p}e^{-i\xi x}V(t,x,k,p).$ (26)

Using Eq. (25), it is simple to see that $\tilde{V}$ obeys

$\displaystyle
i\partial_{t}\tilde{V}=(k\cdot\xi-\eta\Sigma(k))\tilde{V}+\eta\int
dk^{\prime}\sigma(k,k^{\prime})\tilde{V}(t,\xi,k^{\prime},\eta).$ (27)

We discretise in $\xi,k$ and $\eta$ with $\Delta\xi=2/J$, $\Delta k=2/K$ and
$\Delta\eta=2/N$ for positive even integers $J,K,N$. Let
$\boldsymbol{\tilde{V}}(t)=\boldsymbol{\tilde{V}}_{[\xi]}\otimes\boldsymbol{\tilde{V}}_{[k]}\otimes\boldsymbol{\tilde{V}}_{[p]}$
be the vector formed by the discretisation of $\tilde{V}$ in $\xi,k$ and
$\eta$. Here each element of $\boldsymbol{\tilde{V}}_{[x]}$ denotes
$\tilde{V}(t,\xi_{i},k,\eta)$ with $i=1,...,J^{d}$ labelling the $\xi$ grid
points, each element of $\boldsymbol{\tilde{V}}_{[k]}$ is
$\tilde{V}(t,\xi,k_{j},\eta)$ with $j=1,...,K^{d}$ labelling the $k$ grid
points and the elements of $\boldsymbol{\tilde{V}}_{[\eta]}$ are
$\tilde{V}(t,\xi,k,\eta_{k})$ with $k=1,...,N$. One can then use the
quadrature rule for the last term of Eq. (27) to write

$\displaystyle\int
dk^{\prime}\sigma(k,k^{\prime})\tilde{V}(t,\xi,k^{\prime},\eta)\approx\sum_{j^{\prime}=1}^{K^{d}}\sigma_{j,j^{\prime}}\tilde{V}(t,\xi,k^{\prime}_{j^{\prime}},\eta)$
(28)

where $\sigma_{j,j^{\prime}}=\sigma(k_{j},k^{\prime}_{j^{\prime}})/K^{d}$ are
the elements of a symmetric (and Hermitian) matrix $\boldsymbol{\sigma}$,
since we are looking at isotropic scattering. Then from Eq. (27) the vector
$\boldsymbol{\tilde{V}}(t)$ obeys

$\displaystyle
i\frac{d\boldsymbol{\tilde{V}}}{dt}=\boldsymbol{H}_{\text{total}}\boldsymbol{\tilde{V}}$
(29)

where the matrix $\boldsymbol{H}$ is Hermitian

$\displaystyle\boldsymbol{H}_{\text{total}}=\boldsymbol{L}_{[\xi,k]}\otimes\boldsymbol{1}-\mathbf{1}\otimes\boldsymbol{\Sigma}\otimes\boldsymbol{D}+\mathbf{1}\otimes\boldsymbol{\sigma}\otimes\boldsymbol{D}.$
(30)

Here $\boldsymbol{D}=\text{diag}(\mu_{1},...,\mu_{N})$ with
$\mu_{k}=\pi(k-N/2)$ and
$\boldsymbol{\Sigma}=\text{diag}(\Sigma(k_{1}),...,\Sigma(k_{K^{d}}))$. The
matrix $\boldsymbol{L}_{[\xi,k]}$ acts on
$\boldsymbol{\tilde{V}}_{[\xi]}\otimes\boldsymbol{\tilde{V}}_{[k]}$ and
multiplies its $(i,j)^{\text{th}}$ component it by factor $\xi_{i}\cdot
k_{j}$.

Since $\boldsymbol{H}_{\text{max}}$ is Hermitian, one can use Hamiltonian
simulation methods in Eq. (29) where
$\boldsymbol{\tilde{V}}(t)=\exp(-i\boldsymbol{H}_{\text{total}}t)\boldsymbol{\tilde{V}}(0)$
and one can prepare the quantum state
$|\tilde{V}(t)\rangle=(1/\|\boldsymbol{\tilde{V}}(t)\|)\sum_{i=1}^{J^{d}}\sum_{j=1}^{K^{d}}\sum_{k=1}^{N}\boldsymbol{\tilde{V}}(t)|i,j,k\rangle$
by just evolving the Hamiltonian $\boldsymbol{H}_{\text{total}}$, when given
the initial state $|V(0)\rangle$. Note that
$\|\boldsymbol{\tilde{V}}(t)\|=\|\boldsymbol{\tilde{V}}(0)\|$ since
$\boldsymbol{\tilde{V}}(t)$ evolves by a unitary transformation. This can be
used to prepare the quantum state whose amplitude is proportional to the
average Wigner function of the quantum state in a random medium in the
following way.

Let $\boldsymbol{W}(t)$ be the vector whose elements are $W(t,x_{i},k_{j})$ at
grid points $(x_{i},k_{j})$. Suppose we want to prepare the quantum state
$|W(t)\rangle=(1/\|\boldsymbol{W}(t)\|)\sum_{i=1}^{M^{d}}\sum_{j=1}^{K^{d}}\boldsymbol{W}(t)|i,j\rangle$.
By definition,
$\boldsymbol{V}(t)=(\mathcal{F}^{-1}_{x}\otimes\mathbf{1}\otimes\mathcal{F}^{-1}_{p})\boldsymbol{\tilde{V}}(t)$,
where $\mathcal{F}^{-1}_{(\cdot)}$ is the inverse Fourier transform with
respect to variable $(\cdot)$ and $\boldsymbol{V}(t)$ is the discretisation of
the warped phase transformation $V(t,x,k,p)$. Then
$|V(t)\rangle=(1/\|\boldsymbol{V}(t)\|)\sum_{i=1}^{M^{d}}\sum_{j=1}^{K^{d}}\sum_{k=1}^{N}\boldsymbol{W}(t)\exp(-p_{k})|i,j,k\rangle$.
Fourier transforms do not change norms, so
$\|\boldsymbol{V}(t)\|=\|\boldsymbol{\tilde{V}}(t)\|=\|\boldsymbol{\tilde{V}}(0)\|=\|\boldsymbol{V}(0)\|$.
We can use the projection $\hat{P}$ in Appendix B and show $|W(t)\rangle$ and
$|\tilde{V}(t)\rangle$ are related by

$\displaystyle\hat{P}(\mathcal{F}^{-1}_{x}\otimes\mathbf{1}\otimes\mathcal{F}^{-1}_{p})|\tilde{V}(t)\rangle=\frac{\|\boldsymbol{W}(t)\|}{\|\boldsymbol{V}(0)\|}|W(t)\rangle\sum_{k=1}^{N}e^{-p_{k}}|k\rangle$
(31)
$\displaystyle=\frac{\|\boldsymbol{W}(t)\|\|\exp(-p)\|}{\|\boldsymbol{V}(0)\|}|W(t)\rangle|\exp(-p)\rangle$

where $|\exp(-p)\rangle$ is the state with normalisation constant
$\|\exp(-p)\|\sim\mathcal{O}(\sqrt{N})$, as defined in Appendix B. Then
following similar methods in Appendix B, it can be proved that
$\|\boldsymbol{V}(0)\|/\|\boldsymbol{W}(t)\|\sim\sqrt{N}\|\boldsymbol{W}(0)\|/\|\boldsymbol{W}(t)\|$.

The proof of the resource cost in obtaining $|W(t)\rangle$ follows in the same
way as shown in Appendix B. Let $s$ be the sparsity of
$\boldsymbol{H}_{\text{max}}$ and
$\|\boldsymbol{H}_{\text{total}}\|_{\text{max}}$ is its max-norm. Then
following Lemma 2, we see that, given $|W(0)\rangle$, the query and gate
complexities to obtain $|W(t)\rangle$ are respectively
$\tilde{\mathcal{O}}((\|\boldsymbol{W}(0)\|/\|\boldsymbol{W}(t)\|)st\|\boldsymbol{H}_{\text{total}}\|_{\text{max}})$
and
$\tilde{\mathcal{O}}((\|\boldsymbol{W}(0)\|/\|\boldsymbol{W}(t)\|)dst\|\boldsymbol{H}_{\text{total}}\|_{\text{max}})$.

From Eq. (30), one sees that the sparsity $s$ is dominated by
$s(\boldsymbol{\sigma})$, which is the sparsity of $\boldsymbol{\sigma}$. This
means the more off-diagonal terms in $\boldsymbol{\sigma}$ (corresponding to
more scattering to different wave vectors), the larger the cost in preparing
$|W(t)\rangle$. To find the max-norm
$\|\boldsymbol{H}_{\text{total}}\|_{\text{max}}$, we note that
$\|\boldsymbol{L}_{[\xi,k]}\|_{\text{max}}\sim\|\boldsymbol{\sigma}\otimes\boldsymbol{D}\|_{\text{max}}\sim\|\boldsymbol{\sigma}\otimes\boldsymbol{D}\|_{\text{max}}\sim\mathcal{O}(1/\epsilon)$
where $\boldsymbol{\sigma}$ and $\boldsymbol{\Sigma}$ are order one factors.
This means that the Schrödingerisation approach does not change the order of
$\|\boldsymbol{H}_{\text{total}}\|_{\text{max}}$, since the transformed
equation remains first order as the original equation.

We can also apply this to finding the stationary state of the semiclassical
limit of the Wigner function, i.e., $|W(t_{s})\rangle$ when
$d\boldsymbol{W}(t)/dt|_{t=t_{s}}=0$. The analysis proceeds in the same way to
preparing the ground state in Appendix C. The transport operator in Eq. (23)
is dissipative, or rather hypocoercivity, with spectral gap shown in [27]. In
the case of isotropic scattering, the spectra includes $0$ and all other
strictly negative ones, depending on Chandrasekhar’s H-function [28]. We can
thus use $E_{0}=0$ in this case and the same cost as in Appendix C follows.

We can also extend the above analysis in the high $d$ limit by reinterpreting
as a multiparticle $n=d/3$ limit with each particle moving in three spatial
dimensions. However, in this case, $\omega^{2}=k^{2}/2+U_{C}$ since the
interparticle interactions with a potential $U_{C}$ (e.g. Coulomb) needs to be
included. This means we need to augment Eq. (24) by an extra term

$\displaystyle\partial_{t}W+k\cdot\nabla_{x}W-\nabla_{x}U_{C}\cdot\nabla_{k}W$
$\displaystyle=\int
dk^{\prime}\sigma(k,k^{\prime})W(t,x,k^{\prime})-\Sigma(k)W.$ (32)

Here we omit the precise details (see [24]), but one can discretise in $x$,
$k$ and rewrite this as a system of linear ODEs and apply the general
formalism.


In industrial model predictive controllers (MPCs), models generated from regression-based system identification methods typically contain small or even physically non-existent degrees of freedom.  Control issues can arise when the steady-state optimizer uses these small degrees of freedom to calculate targets for plant operation due to matrix ill-conditioning. Mathematical techniques like Relative Gain Array (RGA) and Singular Value Decomposition (SVD) are helpful for analyzing controller gain interactions and identifying conditioning issues, which can be corrected relatively easily in small models.  However, these techniques are difficult and tedious to apply for larger, more complex models. This paper describes a novel, non-iterative, RGA-based, binning technique for discretizing the gain matrix and quickly solving 2$\times$2 conditioning issues for any model size, while guaranteeing gain adjustments below a certain threshold. Higher order interactions are also discussed.

Process Control & Automation, Controller Performance Evaluation, Identification & Estimation

§ INTRODUCTION

Industrial model predictive controllers (MPCs) are commonly implemented with an integrated steady-state, Linear Program (LP) optimizer to generate targets that drive the process to an economic optimum <cit.>. In industrial MPC projects, one of the most difficult and time-consuming steps is understanding the interaction of the steady-state model gains with the optimizer costs in the LP solution. Part of this consideration is respecting the material and energy balance of the process, and determining how much “detail" in terms of smaller gains should be built into the model.

In practice, designers must make structural decisions in the gain matrix to reconcile the model with the true process. This step relies on the designer's process knowledge and experience, along with the usage of process data to justify any curve or gain manipulation of the model. Adjusting the gain matrix is often a balancing act between maintaining model accuracy and usability. Deleting a significant gain or multiplying a gain by an arbitrary factor to try to “tune" the controller should be avoided. A case study from Mitsubishi demonstrated how the erroneous removal of a gain element, thought to be insignificant based on engineering judgment, resulted in controller instabilities <cit.>.

Fundamentally, there are 3 types of MV-CV (Manipulated Variable-Controlled Variable) interactions:

* Non-collinearities: MV-CV interactions that are easy to independently control

* Collinearities: MV-CV interactions that cannot be independently controlled, because the gains are collinear

* Near-collinearities: MV-CV interactions that are difficult to independently control

The non-collinear MV-CV interactions are typically trivial to model and tune, and do not warrant in-depth discussions for the purposes of this paper.

Process flow diagram (PFD) for a debutanizer column, which separates butane and light-end components from a heavier (C5+) hydrocarbon mixture. MVs are highlighted in blue, while CVs are red.

The collinear MV-CV interactions can often be identified based on process data and a priori process knowledge. These interactions may not be exactly collinear in the raw model, but can be corrected by forcing the gain ratios of the collinear MV-CV pairs to be numerically equal. Doing so ensures that the LP will not try to use control handles that are nonexistent.

The near-collinear MV-CV pairs can result in problematic interactions due to numerically ill-conditioned matrices, which, if active, cause undesirable large MV movements to control the near-collinear CVs independently at their constraint limits. In industry, it is common for control practitioners to adjust the steady-state gains manually <cit.> to deal with ill-conditioning issues. When a gain is changed, its interaction with all of the other gains also changes. Consequently, fixing one gain conditioning issue may create other issues in other parts of the matrix. For further details, Waller and Waller (1995) <cit.> provides an excellent review of process directionality literature.

An ExxonMobil patent <cit.> describes a mathematical procedure for adjusting gain matrices. Some commercial applications try to automatically correct the ill-conditioning, such as Honeywell's RMPCT which uses Singular Value Thresholding (SVT) <cit.>, but automatic techniques can fail because the raw gains can be far enough from the collinearity threshold that they are not fixed. On the other hand, making the SVT threshold large enough to fix these gains could incorrectly remove small but critical interactions/degrees of freedom. Semi-automated techniques, such as those described in a patent by Aspentech <cit.>, may not always produce satisfactory results that are consistent with the process, and may warrant further manual adjustments.

In this paper, we present a novel, non-iterative, binning approach to gain conditioning that significantly reduces the difficulties associated with gain matrix adjustments faced by control practitioners.

§ REVIEW OF TECHNIQUES

Two primary tools for analyzing model conditioning are Singular Value Decomposition (SVD), and the Relative Gain Array (RGA).

§.§ Singular Value Decomposition

Consider a $2\times2$ steady-state gain matrix, $G$,

\begin{equation}
\label{eqn:ssgain-general}
G =
\begin{bNiceMatrix}[first-row, last-col=3]
{\footnotesize{\textrm{MV1}}} &  {\footnotesize{\textrm{MV2}}} & \\
G_{11} & G_{12} & \footnotesize{\textrm{CV1}}  \\
G_{21} & G_{22} & \footnotesize{\textrm{CV2}}
\end{bNiceMatrix}
\end{equation}

SVD decomposes $G$ into three matrices with $G = U \Sigma V^\mathsf{T}$, describing its most responsive and least responsive directions in terms of both MV moves and CV responses, where

\begin{equation*}
\Sigma =
\begin{bNiceMatrix}
\overline{\sigma} & 0  \\
0 & \underline{\sigma}
\end{bNiceMatrix} \quad
U =
\begin{bNiceMatrix}[first-row]
{\overline{\bm{u}}} &  {\underline{\bm{u}}} \\
u_{11} & u_{12} \\
u_{21} & u_{22}
\end{bNiceMatrix}
\quad
V^\mathsf{T} =
\begin{bNiceMatrix}[last-col]
v_{11}^\mathsf{T} & v_{12}^\mathsf{T} & \overline{\bm{v}}^\mathsf{T} \\
v_{21}^\mathsf{T} & v_{22}^\mathsf{T} & \underline{\bm{v}}^\mathsf{T}
\end{bNiceMatrix}
\end{equation*}

A physical interpretation of SVD in the context of multivariable control is described in detail by Moore (1986) <cit.>, as illustrated in Fig. <ref>. CVs will be most responsive in the direction of $\overline{\bm{u}}$ when the MVs are moved in the direction of $\overline{\bm{v}}^\mathsf{T}$. The opposite is true for $\underline{\bm{u}}$ and $\underline{\bm{v}}^\mathsf{T}$. The singular values $(\overline{\sigma}, \underline{\sigma})$ represent the magnitude of each CV direction's controllability.

Singular values are proportional to system sensitivity <cit.>; high $\sigma$ values indicate that MV moves will need to be small to maintain accurate control, and vice versa. Singular values that are close to each other indicate that, generally, a controller will be able to control different CVs relatively independently. In contrast, systems with very different singular values will require larger MV moves for disproportionately smaller CV changes. This effect is encapsulated by the condition number, which is defined as the ratio of the largest to the smallest singular value:

\begin{equation}
\gamma = {\overline{\sigma}}/{\underline{\sigma}}
\end{equation}

It is important to note that the SVD is scale-dependent. The singular values and condition numbers calculated are based on the engineering units of the MVs and CVs.

As an example, consider an unscaled gain matrix in Eq. (<ref>) based on the debutanizer process in Figure <ref>. MV1 is the , MV2 is the , CV1 is the  in the bottoms and CV2 is the  in the overhead.

\begin{equation}
\label{eqn:ssgain}
G =
\begin{bNiceMatrix}[first-row, last-col=3]
{\footnotesize{\textrm{MV1}}} &  {\footnotesize{\textrm{MV2}}} & \\
-0.1942 & -0.0029 & \footnotesize{\textrm{CV1}}  \\
0.1843 & -0.0288 & \footnotesize{\textrm{CV2}}
\end{bNiceMatrix}
\end{equation}

The SVD results based on Eq. (<ref>) are summarized below.

\begin{equation*}
U =
\begin{bmatrix*}[r]
-0.7215 & 0.6925 \\
0.6925 & 0.7215
\end{bmatrix*}
\quad
V =
\begin{bmatrix*}[r]
0.9978 & -0.0664 \\
-0.0664 & -0.9978
\end{bmatrix*}
\end{equation*}

\begin{equation*}
\Sigma =
\begin{bmatrix*}[c]
0.2683 & 0 \\
0 & 0.0229
\end{bmatrix*}
\end{equation*}

The condition number is calculated to be $\gamma = 11.74$. However, due to the scale-dependent properties of SVD, the results obtained with Eq. (<ref>) may not be meaningful.

Graphical interpretation of SVD. The left subfigure represents a unit circle of MV input moves, and the right subfigure represents the CV responses. SVD breaks down the CV responses in terms of a `strong' (red) and `weak' (blue) direction illustrated by the CV ellipse. As the weak direction becomes weaker, $\gamma$ increases and the system converges towards collinearity. A collinear matrix would collapse the CV ellipse into a line in the direction of $\overline{\bm{u}}$.

Researchers have proposed various scaling choices, including scaling by physical limits, relative importance of variables, or by other problem-specific scaling that minimizes the condition number to obtain $\gamma_{\min}$ <cit.>. In practice, a common and useful way to scale the matrix is typical move scaling.
This method multiplies each gain by the maximum MV move size made during identification, then divides it by the maximum gain of that CV across all other MVs. The scaled matrix, $G'$ using typical move scaling is given by

\begin{equation}
% Notation used in Waller and Waller
% https://math.stackexchange.com/questions/3274790/notation-for-element-wise-multiplication-of-vector-and-matrix-columns
\label{eqn:svd_scaling}
G' = S_b G S_a
\end{equation}

where $S_a$ and $S_b$ are non-zero, diagonal scaling matrices such that
S_a = \diag\left(\Delta\textrm{MV}_1, \Delta\textrm{MV}_2\right)
S_b = \diag(\underset{i}{\max} (|GS_{a}|)_{1, i}, \underset{i}{\max} (|GS_{a}|)_{2, i})^{-1}
In the first scaling step, the diagonal elements in $S_a$ scale each MV column in $G$ by its $\Delta\textrm{MV}$ move size. In the second scaling step, the diagonal elements in $S_b$ scale each CV row by dividing each row by that CV's maximum gain in the $GS_a$ matrix.

Typical move scaling `normalizes' the gains by considering both MV move sizes and CV responses to all related MVs. Consequently, small scaled gains may indicate that the gain is small or insignificant, or the MV move size was too small for that CV, or that large MV moves were made with a high-gain CV, making other CV scaled gains look small.

We now consider a scaled gain matrix based on Eq. (<ref>). Let the  move size, $\Delta\textrm{MV}_1 = 2$, and the  move size, $\Delta\textrm{MV}_2 = 10$. Using typical move scaling, the gains in each CV now lie between -1 to +1. The scaled gain matrix, $G'$ is calculated using Eq. (<ref>) as

\begin{equation}
\label{eqn:ssgain_scaled}
G' =
\begin{bNiceMatrix}[first-row, last-col=3]
{\footnotesize{\textrm{MV1}}} &  {\footnotesize{\textrm{MV2}}} & \\
-1.0000 & -0.0754 & \footnotesize{\textrm{CV1}}  \\
1.0000 & -0.7813 & \footnotesize{\textrm{CV2}}
\end{bNiceMatrix}
\end{equation}

Scaling has changed our understanding of the problem. The condition number, $\gamma' = 2.68$ is now smaller for $G'$, compared to $\gamma = 11.7$ for $G$ prior to scaling. The higher the condition number, the harder it is to drive in the weak direction. From a distillation control perspective, we know that it is much easier to lower the RVP and increase the C5s simultaneously (shift the column composition profile u p), than it is to lower both the RVP and C5s (increase separation) <cit.>.

§.§ Relative Gain Array

Inspecting the RGA value is another approach to identify gain ratio issues <cit.>.  The advantage of RGA over SVD is that it is scale independent. In practice, for the gains in a multivariable controller matrix, the RGA is useful in setting an interaction threshold for significant and insignificant gain ratios. For the case of the $2\times2$ matrix, $G$ in Eq. (<ref>), the RGA computations can be simplified <cit.> as:

\begin{equation}
\label{eqn:RGA}
\textsc{RGA} = \left(1 - \frac{G_{12}\cdot G_{21}}{G_{11}\cdot G_{22}}\right)^{-1}
\end{equation}

As the RGA number increases, the gain ratios in the $2\times2$ matrix approach each other and the matrix approaches singularity.  When the gains are collinear, the RGA number is infinite. For an APC model, it is useful to set an RGA threshold such that any gain ratios above a threshold RGA value need to be fixed.  For $G$ in Eq. (<ref>), the RGA is calculated to be a desirable, small value of 1.11. In industrial practice, a reasonable 2x2 RGA threshold in a typical multivariable controller could be around 12. Because RGA is independent of scaling, the RGA describes process directionality with regard to inherent process stability, while SVD, because it is considering gain magnitude, does so with regard to performance <cit.>. From an MPC standpoint, we may or may not be concerned from the performance point of view. Economically we still may want to use the MV with a small directional effect because it costs nothing to do so. For further details and control applications of RGA, the expository text by McAvoy <cit.> is an excellent resource.

§ METHODOLOGY: BINNING THE GAINS

We now present a novel, non-iterative methodology for gain matrix conditioning. The traditional approach involves tedious manual gain adjustments or a complicated optimization analysis that may or may not converge, until a reasonable RGA value is reached for all submatrices. Our methodology takes a user-defined RGA threshold to generate a discretized grid of gain values with desirable properties for matrix conditioning.

Consider any $2\times2$ non-zero gain matrix, $G$. The gains can be scaled to obtain $K$, such that 3 of the 4 gains are $1$, with a single non-unity gain, $k$, that is dependent on the matrix's RGA value:

\begin{equation}
\label{eqn:unityscaling}
\begin{aligned}
K &=
\begin{bmatrix*}[c]
1/G_{11} & 0 \\
0 & 1/G_{21}
\end{bmatrix*}
\begin{bmatrix*}[c]
G_{11} & G_{12} \\
G_{21} & G_{22}
\end{bmatrix*}
\begin{bmatrix*}[c]
1 & 0 \\
0 & G_{21}/G_{22}
\end{bmatrix*} \\
& = \begin{bmatrix*}[c]
1 & k \\
1 & 1
\end{bmatrix*}, \quad \text{where } k = (1-\RGA^{-1})
\end{aligned}
\end{equation}

The $k$ term defines the relationship between the gain ratios and the RGA, as shown by rearranging Eq. (<ref>) to get

\begin{eqnarray}
G'_{12}/G'_{22} = (1 - \textsc{RGA}^{-1}) \cdot G'_{11}/G'_{21} \label{eqn:RGA2_col}\\
G'_{12}/G'_{11} = (1 - \textsc{RGA}^{-1}) \cdot G'_{22}/G'_{21} \label{eqn:RGA2_row}
\end{eqnarray}

where Eq. (<ref>) defines the relationship for the column ratios, and likewise, Eq. (<ref>) for the rows. For $\RGA = \infty$, $k=1$ and the matrix is singular, meaning that physically, either the MVs or CVs are collinear. For $\RGA = 12$, $k = 0.9167$. The intuition behind our method is that, instead of just using the matrix gains to compute $k$ and $\RGA$, we could also flip the relationship around to enforce a certain $\RGA$ value for any matrix. This is done by defining our desired $\RGA$, calculating $k$ and then adjusting the gains to equal $k$.

In any typical-move scaled matrix obtained with Eq. (<ref>), the absolute values of the scaled gains lie between $0$ and $1$. Starting with the largest scaled gain of $1$, for $\RGA=12$, the next-largest value is $k=0.9167$ from Eq. (<ref>). To `condition' the matrix, gains between $1$ and $0.9167$ are adjusted to the closest value, either to $1$, enforcing collinearity, or to $0.9167$, enforcing $\RGA = 12$.

To cover the entire range of typical move-scaled gains, we now seek a series of values, $k_1, \cdots, k_n$ such that all scaled gains can be adjusted to those values, and the entire matrix can be conditioned to any user-defined RGA threshold. We now ask, what is the next largest value of $k$ that will ensure that the matrix meets the desired $\RGA=12$ threshold? It is found by multiplying $K$ by $k$:

\begin{equation}
\label{eqn:unityscaling2}
\begin{aligned}
K^{(2)} &= k\cdot \begin{bmatrix*}[c]
1 & k \\
1 & 1
\end{bmatrix*}
= \begin{bmatrix*}[c]
0.9167 & 0.8403 \\
0.9167 & 0.9167
\end{bmatrix*}
\end{aligned}
\end{equation}

Similar to the first step, the absolute value of the scaled gains between $0.9167$ and $0.8403$ are adjusted to the closest value, either $0.9167$ or $0.8403$. We can repeat this procedure to get $K^{(3)}, \cdots, K^{(n)}$ until we generate enough values to cover the range of the scaled gain matrix.

Formally, the pattern in Eq. (<ref>) can be used to construct a recurrence relation that yields a geometric sequence of gains that satisfy a desired RGA threshold. This sequence forms the basis of a set of discretized bins. We define bins as having strictly decreasing bin boundaries, $B_0, B_1, \cdots, B_n$ such that $B_{i+1} < B_{i}$ for $i = {0,1,\cdots,n}$. The bin boundaries define the gain values we can use. The bin width for each interval, $[B_{i+1}, B_{i}]$, is given by $w_i = B_{i} - B_{i+1}$. This definition constructs $n$ bins with $n+1$ bin boundaries. $n$ is chosen such that the bins cover the range of all gains in the matrix.

Let the first bin boundary be $B_{0} = 1$. The next term, $B_1$, is found with a recurrence relation derived from Eq. (<ref>) as

$$B_1 = (1 - \textsc{RGA}^{-1}) \cdot B_{0}$$

Solving the recurrence relation yields a closed-form solution for any $i$th bin boundary as a function of $B_0$ and RGA:

\begin{equation}
\label{eqn:RGA_buckets}
B_i = (1 - \textsc{RGA}^{-1})^i \cdot B_0
\end{equation}

Consider a scaled gain element $G'_{ij}$ in $G'$, with an absolute value $|G'_{ij}|$ that lies in the interval $B_{k+1} \leq |G'_{ij}| < B_{k}$. The binning procedure discretizes $G'_{ij}$ to obtain a binned gain, $\widehat{G}_{ij}$, by assigning $G'_{ij}$ to the closest boundary with

\begin{equation}
\label{eqn:binning_step}
\widehat{G}_{ij}=
\sgn{(G'_{ij})} \times
\begin{cases}
B_{k},& \text{if } |G'_{ij}| > \frac{B_k + B_{k+1}}{2}\\
B_{k+1},              & \text{otherwise}
\end{cases}
\end{equation}

In other words, the procedure simply adjusts the typical-move scaled gain to the closest bin boundary, based on the absolute value of the gain. Due to the exponential term in Eq. (<ref>), the bin widths are also strictly decreasing, such that $w_{k+1} < w_{k}$. It follows that the maximum gain change introduced by this procedure is bounded by the width of the first bin interval $[B_1, B_0]$, with the largest gain adjustment occurring at the midpoint of $[B_1, B_0]$, denoted by $B_m$. As shown by the magenta marker in Figure <ref>, a gain at the midpoint of $[B_6, B_5]$ could be assigned to either boundaries and receive an equal amount of gain adjustment.

axis line style=thick, cyan,
axis x line=middle,
xticklabel style = yshift=0cm,
every x tick/.style=thick, color=teal,
major tick length=0.20cm,
axis y line=none,
[only marks,mark=*,red] coordinates
(0.7813,0) (1.0,0)
[only marks,mark=*,red] coordinates
[label=180:Scaled, $g'$](x0) at (axis cs:0.48,0);
[pin=[pin edge=thin, pin distance=0.4cm] above right:1.000](x1) at (axis cs:1,0);
[pin=[pin edge=thin, pin distance=0.4cm] above right:0.7813](x2) at (axis cs:0.7813,0);
[pin=[pin edge=thin, pin distance=0.4cm] above left:0.6203](x3) at (axis cs:0.6203,0);

ticks = none,
axis line style=thick, cyan,
axis x line=middle,
major tick length=0.25cm,
axis y line=none,
[only marks,mark=|,thick,black,mark size=0.125cm] coordinates
(0.5439,0) (0.5933,0) (0.6472,0) (0.7061,0) (0.7703,0) (0.8403,0) (0.9167,0) (1.0,0)

[only marks,mark=diamond*,magenta] coordinates

[label=180:Binned, $\hat{g}$](y0) at (axis cs:0.48,0);

[label=[label distance=]270:$B_7$](b7) at (axis cs:0.5439,0);

[label=[label distance=]270:$B_4$](b7) at (axis cs:0.7061,0);
[label=[label distance=]270:$B_3$](b7) at (axis cs:0.7703,0);
[label=[label distance=]270:$B_2$](b7) at (axis cs:0.8403,0);
[label=[label distance=]270:$B_1$](b7) at (axis cs:0.9167,0);
[label=[label distance=]270:$B_0$](b7) at (axis cs:1,0);

(ymid) at (axis cs:0.6203,0) ;
(ymidL) at (axis cs:0.5933,0) ;
(ymidR) at (axis cs:0.6472,0) ;

[pin=[pin edge=thin, pin distance=0.4cm] above right:1.000](y1) at (axis cs:1,0);
[pin=[pin edge=thin, pin distance=0.4cm] above right:0.7703](y2) at (axis cs:0.7703,0);
[pin=[pin edge=thin, pin distance=0.4cm] above left:0.5933](y4) at (axis cs:0.5933,0);
[pin=[pin edge=thin, pin distance=0.4cm] above right:0.6472](y3) at (axis cs:0.6472,0);

[-latex, dotted, shorten >=0.15cm, semithick] (x1) – (y1);
[-latex, dotted, shorten >=0.15cm, semithick] (x2) – (y2);
[-latex, dotted, magenta, shorten >=0.15cm, semithick] (x3) – (y3);
[-latex, dotted, magenta, shorten >=0.15cm, semithick] (x3) – (y4);

[to-to, magenta] ([yshift=-0.2cm]ymid.center) – ([yshift=-0.2cm]ymidR.center) node[midway,below] $\frac{w_5}{2}$;

[to-to, magenta] ([yshift=-0.2cm]ymid.center) – ([yshift=-0.2cm]ymidL.center);

Illustration of the binning process. A gain, $g'$ that lies in the bin interval $[B_{i+1}, B_i]$ is adjusted to the closest bin boundary to output a binned gain, $\hat{g}$. Gains that lie on the midpoint of an interval are at equidistant from both bin boundaries, and will experience equal gain adjustment in either direction, as shown by the magenta arrows.

The maximum relative gain change (%), $\delta_{\max}$ is a function of only RGA and independent of the choice of $B_0$. By definition, the relative gain change is the ratio of the gain adjustment to the initial gain value, $B_m$

\begin{equation}
\label{eqn:relative_change}
\delta_{\max} = \frac{B_0 - B_m}{B_m}\times 100\% \text{,} \quad B_m = \frac{B_0+B_1}{2}
\end{equation}

By substituting Eq. (<ref>) into Eq. (<ref>), the $B_0$ terms cancel out and we get $
\delta_{\max} = \frac{\RGA^{-1}}{2 - \RGA^{-1}} \times 100\%$. The maximum relative gain change percentage as a function of RGA choice is shown in Figure <ref>. For an RGA choice of 12, the largest relative gain change percentage would have an upper bound of 4.34%, where the 0.9583 midpoint gain is adjusted to either $B_1=0.9167$ or to $B_0 = 1.0$.

xmin=0, xmax=20,
ymin=0, ymax=25,
minor y tick num=4,

xlabel=RGA choice,
ylabel=$\delta_{\max}$ (%)
[only marks,black,samples at=12] 100*abs((1/x)/((1/x)-2)) node[pin=30:$(12,4.35\%)$];
[only marks,black,samples at=6] 100*abs((1/x)/((1/x)-2)) node[pin=30:$(6,9.09\%)$];
[anchor=south west, text=black] (a) at (axis cs:8,17) $\delta_{\max} = \dfrac{\RGA^{-1}}{2 -\RGA^{-1}} \times 100\%$;

The maximum relative gain change increases exponentially with smaller RGAs, showing that larger gain adjustments are required to achieve lower RGA values. In practice, for a reasonable RGA threshold of 12, any gain change is guaranteed to at most 4.35%.

The selected RGA threshold builds a binning structure for the typical move scaled model gains to be fitted into. There is no iteration, and when new gains are added to the gain matrix, they are adjusted to the existing structure, rather than iteratively adjusting the entire matrix of gains. It is also important to note that only the gains that cause $2 \times 2$ matrices to exceed the $\textsc{RGA}$ threshold need to be adjusted, and all others can be left as-is.

The binning procedure is applied to the scaled gain matrix, $G'$ in Eq. (<ref>) to obtain binned gains $\widehat{G}$ as shown in Eq. (<ref>). CV1-MV2 was adjusted by 2.5% and CV2-MV2 by 1.4%.

\begin{equation}
\label{eqn:ssgain_binned}
\widehat{G} =
\begin{bNiceMatrix}[first-row, last-col=3]
{\footnotesize{\textrm{MV1}}} &  {\footnotesize{\textrm{MV2}}} & \\
-1.0000 & -0.07351 & \footnotesize{\textrm{CV1}}  \\
1.0000 & -0.77025 & \footnotesize{\textrm{CV2}}
\end{bNiceMatrix}
\end{equation}

§ CASE STUDY: DEBUTANIZER PROBLEM

We apply these techniques on the debutanizer process described in Figure <ref>. This process is a good pedagogical example — although the operation of a debutanizer seems relatively simple, the process illustrates many complexities that are commonly faced by APC engineers. The raw gain matrix for the process is shown in Table <ref>. Each of these variables is shown in the debutanizer schematic in Fig. <ref>.

Original Gains, $G$

TC-REBOIL-SP     FC-REFLUX-SP     PC-TOP-SP     FC-DIST-SP     FI-FEED-PV

AI-RVP-PV         $-0.1942$          $-0.0029$          $0.0711$           $0$               $0.0013$

AI-DIST-C5         $0.1843$           $-0.0288$          $-0.1907$          $0$               $0.0070$

TOP-PCT           $0.9220$           $-0.1477$          $-0.9458$          $0$               $0.0371$

LI-ACCUM-PF       $0.2042$           $0$                $-0.0667$          $-0.1485$         $0.0381$

DP-DEBUT-PV       $0.0774$           $0.0063$           $-0.0143$          $0$               $0.0064$

PC-TOP-OPT        $4.9714$           $0.5000$           $-4.9887$          $0$               $0.3738$

FC-REBOIL-OP      $4.5005$           $0.3391$           $-1.4486$          $0$               $0.2725$

FC-REFLUX-OP      $0$                $0.2651$           $0$                $0$               $0$

The raw gains are scaled according to their typical moves using Eq. (<ref>), and analyzed using SVD and RGA to identify ill-conditioned $2\times2$ submatrices based on an RGA threshold of 12 and a corresponding condition number of 59, as shown in Table <ref>. The markers indicate that the gain is part of an ill-conditioned $2\times 2$ submatrix based on SVD results (red stars, ), and also on RGA results (blue squares, ). As shown in Table <ref>, there are 13 pairs above the condition number threshold, and 11 pairs above the RGA threshold. Notice that the first entry with a condition number of 59.1 has a corresponding RGA number of 14.4, while another pair with a condition number of 60 has an RGA number of 10.75. For higher-order sub-matrices, there are 34 $3 \times 3$ submatrices of gains with a condition number greater than 100, and 36 $4 \times 4$ submatrices with a condition number greater than 100.

Scaled Gains, $G'$

TC-REBOIL-SP     FC-REFLUX-SP     PC-TOP-SP     FC-DIST-SP     FI-FEED-PV

$\Delta\textrm{MV} = 2$     $\Delta\textrm{MV} = 10$     $\Delta\textrm{MV} = 2$     $\Delta\textrm{MV} = 5$     $\Delta\textrm{MV} = 10$

AI-RVP-PV        $-1$          $-0.0754$     $0.3664$      $0$       $0.0337$

AI-DIST-C5        $0.9666$      $-0.7552$     $-1$     $0$       $0.1839$

TOP-PCT          $0.9748$      $-0.7807$     $-1$     $0$       $0.1962$

LI-ACCUM-PF      $0.5500$      $0$           $-0.1797$     $-1$      $0.5129$

DP-DEBUT-PV      $1$           $0.4049$      $-0.1848$     $0$       $0.4145$

PC-TOP-OPT       $0.9965$      $0.5011$      $-1$     $0$       $0.3747$

FC-REBOIL-OP     $1$           $0.3767$      $-0.3219$     $0$       $0.3027$

FC-REFLUX-OP     $0$           $1$           $0$     $0$       $0$

Ill-conditioned pairs in $G'$

MV1        MV2           CV1        CV2      $\gamma$     RGA

FC-REFLUX-SP     FI-FEED-PV       PC-TOP-OPT       FC-REBOIL-OP     59.14      14.36

TC-REBOIL-SP     FI-FEED-PV       DP-DEBUT-PV      PC-TOP-OPT       59.99      10.75

TC-REBOIL-SP     PC-TOP-SP        AI-RVP-PV        LI-ACCUM-PF      67.50      9.26

TC-REBOIL-SP     FC-REFLUX-SP     DP-DEBUT-PV      FC-REBOIL-OP     81.83      14.37

FC-REFLUX-SP     PC-TOP-SP        AI-DIST-C5        TOP-PCT          124.38     30.54

TC-REBOIL-SP     PC-TOP-SP        AI-DIST-C5        PC-TOP-OPT       131.01     33.24

TC-REBOIL-SP     FC-REFLUX-SP     AI-DIST-C5        TOP-PCT          165.64     40.79

PC-TOP-SP        FI-FEED-PV       AI-DIST-C5        TOP-PCT          169.40     16.04

TC-REBOIL-SP     PC-TOP-SP        TOP-PCT          PC-TOP-OPT       181.27     45.81

TC-REBOIL-SP     FI-FEED-PV       AI-DIST-C5        TOP-PCT          189.76     18.39

FC-REFLUX-SP     FI-FEED-PV       AI-DIST-C5        TOP-PCT          276.03     32.66

TC-REBOIL-SP     PC-TOP-SP        AI-DIST-C5        TOP-PCT          472.37     118.54

TC-REBOIL-SP     PC-TOP-SP        LI-ACCUM-PF      FC-REBOIL-OP     530.00     66.23

Applying the binning strategy of Eq. (<ref>) to the scaled gains yields the binned gains in Table <ref>, which shows each gain's percentage change in red. The resulting gain matrix has three $2 \times 2$ submatrices with an SVD condition number above 59. To eliminate the pairs identified by the condition number threshold, larger gain changes would be required. One of the SVD-identified pairs has an RGA of 8.4. It is therefore unreasonable to remove or adjust the gains associated with this pair because of a somewhat arbitrary SVD condition number. Waller and Waller <cit.> came to the same conclusion related to using the condition number as a measure of ill-conditioning decades ago.

Binned Gains, $\widehat{G}$

TC-REBOIL-SP     FC-REFLUX-SP     PC-TOP-SP     FC-DIST-SP     FI-FEED-PV

AI-RVP-PV        $-1$                       $-0.0754$                       $0.3664$                    $0$           $0.0337$

AI-DIST-C5        $1$ 3.46              $-0.7703$ 2.00             $-1$                        $0$           $0.1914$ 4.08

TOP-PCT          $1$ 2.59              $-0.7703$ -1.34            $-1$                        $0$           $0.1914$ -2.41

LI-ACCUM-PF      $0.5439$ -1.11        $0$                             $-0.1755$ -2.37        $-1$          $0.5129$

DP-DEBUT-PV      $1$                        $0.4189$  3.46             $-0.1848$                   $0$           $0.4189$ 1.06

PC-TOP-OPT       $1$ 0.35              $0.4985$  -0.52            $-1$                        $0$           $0.3840$ 2.49

FC-REBOIL-OP     $1$                        $0.3840$  1.93             $-0.3227$ 0.25         $0$           $0.2958$ -2.30

FC-REFLUX-OP     $0$                        $1$                             $0$                         $0$           $0$

Reviewing the matrix, the binned gain matrix has 10 collinear $2 \times 2$ pairs, where there were none in the raw matrix. Table <ref> highlights the new collinear pairings.

Collinear Pairs in $\widehat{G}$

Pair     MV1        MV2           CV1        CV2

1     TC-REBOIL-SP     FC-REFLUX-SP     AI-C5A5INF2L     TP-TOP-LIN

2     TC-REBOIL-SP     PC-TOP-SP        AI-C5A5INF2L     TP-TOP-LIN

3     TC-REBOIL-SP     FI-FEED-PV       AI-C5A5INF2L     TP-TOP-LIN

4     FC-REFLUX-SP     PC-TOP-SP        AI-C5A5INF2L     TP-TOP-LIN

5     FC-REFLUX-SP     FI-FEED-PV       AI-C5A5INF2L     TP-TOP-LIN

6     PC-TOP-SP        FI-FEED-PV       AI-C5A5INF2L     TP-TOP-LIN

7     FC-REFLUX-SP     FI-FEED-PV       PC-TOP-OPT       FC-REBOIL-OP

8     TC-REBOIL-SP     PC-TOP-SP        LI-ACCUM-PF      FC-REBOIL-OP

9     TC-REBOIL-SP     PC-TOP-SP        AI-C5A5INF2L     PC-TOP-OPT

10     TC-REBOIL-SP     PC-TOP-SP        TP-TOP-LIN       PC-TOP-OPT

Notably, the number of $3\times 3$ gain combinations in the binned gain matrix above a condition number of 100 has been reduced from 34 in the raw matrix to two. The condition numbers of the two $3 \times 3$ matrices are 135 and 105. Similarly, the number of $4 \times 4$ gain combinations with a higher condition number than 100 has been reduced from 36 to 1, where that condition number is 156. By conditioning the $2 \times 2$ gains, the conditioning of the entire gain matrix has improved significantly. It is important to note that a high condition number on a higher-order matrix does not necessarily mean that this combination of gains will be active in an LP solution. Still, reducing the high condition number interactions means that it is impossible for such combinations to play a part in the LP solution under any situation.

§ DISCUSSIONS

How do we know if a gain pair should be collinear or not collinear? Often the answer relies on process knowledge. Second to process understanding, if a gain pair is nearly collinear, removing the degree of freedom based on a high RGA or large condition number is reasonable.

In the debutanizer example, the  is a pressure compensated temperature, and it should be collinear with the  composition. In the binned matrix, it is collinear since the scaled gains are equal. However, if one or more of the gains had been slightly bigger or smaller, the binning process would put the gain in a bin above or below the correct bin. If this were to happen, a degree of freedom would exist in the model where none should be. The opposite situation is also possible, where a degree of freedom is removed in a set of gains where one should exist.

Under certain circumstances, ill-conditioned gains may not be a problem in the controller. Why would we even include the  in the model if we have  which serves the same purpose of measuring C5 composition? One reason is the C5 analyzer may not be reliable. If the controller is intermittently hitting the  or  upper limit, with an unconditioned gain matrix, the LP may believe that temporarily it can control both CVs. As long as the CV limits are spread widely, and the gains are nearly collinear, there likely will be no observed problems. The RGA number is high, but the LP is not trading off the degrees of freedom that exist between the two CVs.

Conditioning problems for this case become apparent when the analyzer starts to drift high or low, and the controller reaches a point when it is trying to control one CV at a high limit, and the other at a low limit.  A properly conditioned model will have no difference in the gain ratios of the two CVs, and when the LP is faced with conflicting limits, it will simply give up on the least important CV.  An unconditioned model will find small degrees of freedom that will allow the LP to trade off MV moves that try to fix the problem.  If the CV gain ratios are a little farther from collinearity, the moves to control CVs at different limits may not be large, but the LP will drive the process in unreasonable directions and CV errors will not be reduced. Also, when differences in gain ratios of collinear variables exist, the controller LP may find a different solution based on which CV limit is active.  This can lead to optimization cost tuning in order to get the controller to behave the same when faced with high or low limits on either variable.

The gain issues described can spread to higher dimensions as well, and the issues are often magnified when the optimizer becomes infeasible.  As long as the MV directions that try to solve the infeasibility are the same as the feasible economic direction, little changes may be noticed when the solution becomes infeasible.  When the infeasible direction is in the opposite direction of the economic direction, the controller can make large moves in an uneconomic direction trying to reduce or eliminate CV error.  From these behaviors, the practitioner may erroneously conclude that the model is inaccurate and needs re-evaluation.  A “good" but unconditioned model can account for this behavior. Retesting to fix an ill-conditioned matrix is unlikely to solve the problem.

These simple examples show why binning the gains is not the final answer, rather the first step in the overall model gain analysis. The “correct" structure of the gain matrix must be determined based on plant test data and engineering fundamentals or process understanding. Binning the gains makes the next step of validating the model gain structure much easier. In a real-world MPC project, the practitioner must perform an interaction analysis of the model gains by simulating the LP solution.

§ CONCLUSIONS

In this paper, we have demonstrated a binning procedure to condition the gain matrix used in industrial MPCs. By conditioning the matrix, we remove small degrees of freedom that are either insignificant or zero, without large perturbations to the process gains. This gain conditioning step prevents the optimizer from trying to use small degrees of freedom to control and optimize the process.  The binning approach also allows for the inclusion of small gains in the model to maintain gain consistency and respect material balances, without introducing further conditioning problems.

There are several areas of interesting future work here. The tables presented earlier for identifying gain interactions using markers are quite basic - as these tables grow in size, this colour-based strategy can result in cluttered and unintuitive displays, making near-collinearity analysis visually difficult. There is an opportunity for improving the visualization of degrees of freedom and gain ratios. We have made some recent progress in visualizing the LP solution for troubleshooting and understanding online MPCs <cit.>.

Unlike SVD, higher-order RGA calculations rely on the matrix inverse, which does not exist for a rank-deficient matrix with collinearities. Therefore, to identify conditioning issues in larger submatrices, the RGA is mathematically undefined and cannot be used, and SVD may not be entirely reliable either due to scaling issues. In industry, higher-dimensional analysis of the gain matrix beyond $3\times3$ submatrices are often not investigated by practitioners. It may be that a higher order RGA analysis could be useful once all $2\times2$ submatrices have been conditioned to the desired RGA threshold. The goal, in practice, has always been to certify that a model is conditioned to the dimensionality of the problem.  Here the dimensionality is defined as a square case, corresponding to the lesser of the number of MVs or the number CVs. For a non-square SVD case, the extra variables generally guarantee that there are significant directions related to the degrees of freedom in the matrix. With SVD, it seems reasonable that for a specific problem, we should be able to identify physically meaningful orthogonal directions in the model, and then force the model gains to follow combinations of these directions.