# Computational Analysis of
Control Systems
Using Dynamic Optimization ††thanks: This work was supported in part by U.S.
Naval Research Laboratory and Air Force Office of Scientific Research

Wei Kang and Liang Xu Wei Kang is with Faculty of Applied Mathematics, Naval
Postgraduate School, Monterey, CA, USA<EMAIL_ADDRESS>Xu is at Naval
Research Laboratory, Monterey, CA, USA<EMAIL_ADDRESS>
###### Abstract

Several concepts on the measure of observability, reachability, and robustness
are defined and illustrated for both linear and nonlinear control systems.
Defined by using computational dynamic optimization, these concepts are
applicable to a wide spectrum of problems. Some questions addressed include
the observability based on user-information, the determination of strong
observability vs. weak observability, partial observability of complex
systems, the computation of $L^{2}$-gain for nonlinear control systems, and
the measure of reachability in the presence of state constraints. Examples on
dynamic systems defined by both ordinary and partial differential equations
are shown.

##  1 Introduction

Control systems are analyzed and characterized by using fundamental concepts
such as observability, reachability, and input-to-output gain [1, 2, 3]. These
concepts have a vast volume of literature. For nonlinear systems, the
challenge is to define the concepts so that they are characteristic and
fundamental to control systems and, meanwhile, they are practically
verifiable. In this paper, the goal is to use dynamic optimization to define
quantitative measures of control system properties. Moreover, computational
methods of dynamic optimization provide practical tools to numerically
implement these concepts in applications.

In Section 2, the ambiguity in estimation is defined as a measure of
observability. This quantity can be numerically computed by solving a dynamic
optimization. An example is shown in which some systems observable in
traditional sense are not practically observable because of their poor value
of ambiguity in estimation. In other words, we can quantitatively tell
strongly observable from weakly observable. Another feature of this concept is
the capability of taking into account non-sensor information or user knowledge
of systems, in addition to the output. For instance, an example is shown in
which the system is unobservable under a traditional definition. It turns out
that the system is strongly observable if we know the control input has a
bounded variation, but without using its accurate upper bound. Moreover, this
concept can be used to measure partial observability of complex systems,
including the observability of a function of the states and the observability
of unknown parameters in a model.

In Section 3, computational methods for the $L^{p}$-gain of control systems
are introduced. The assumption is that the space of input has finite
dimension. Then, the $L^{p}$-gain can be computed using dynamic optimization.
In addition, a method of approximating $L^{p}$-gain is also introduced, which
is based on the eigenvalues of covariance matrices. The methods are
exemplified by a nonlinear model of atomic force microscope.

In Section 4, we define the concepts of ambiguity in control and control cost.
These definitions take into account the control input as well as systems’
constraints. For instance, the concept can be used to quantitatively measure
the reachability of nonlinear systems under the constraint that the states
must stay in a given region of safety. As an example, heat equation with
boundary control is studied.

##  2 Observability

Consider a general control system

$\begin{array}[]{lllllllll}\dot{x}=f(t,x,u,\mu),&x\in\Re^{n_{x}},&u\in\Re^{n_{u}}&\mu\in\Re^{n_{\mu}}\\\
y=h(t,x,u,\mu),&y\in\Re^{n_{y}}\\\ z=e(t,x,u,\mu),&z\in\Re^{n_{z}}\\\
(x(\cdot),u(\cdot),\mu)\in{\cal C}\end{array}$ (1)

in which $y$ is the output, $z$ is the variable to be estimated, which is
either the state $x$ or a function of $x$ in the case of partial observability
for large scale systems. The system state is $x$, $u$ is the control input,
$\mu$ is the parameter or model uncertainty. In (1), ${\cal C}$ is a general
formulation of constraints. Some examples of constraints include, but not
limited to,

$\begin{array}[]{lllllllll}E(x(t_{0}),x(t_{f}))\leq 0,&\mbox{ end point
condition }\\\ s(x,u)\leq 0,&\mbox{ state-control constraints }\\\
s(x(t_{1}))=0,&\mbox{ known event at time }t_{1}\\\
\mu_{min}\leq\mu\leq\mu_{max}&\mbox{ model uncertainties }\\\
s(x,\mu)=0,&\mbox{ DAE (differential-algebraic equations, $\mu$ is a variable)
}\\\ \mbox{Variation}(u)\leq V_{max},&\begin{array}[]{ll}\mbox{control input
with bounded variation}\\\ \mbox{(non-sensor
information)}\end{array}\end{array}$

These constraints represent known information about the system in addition to
the measured output $y$. This general form of constraints makes it possible to
take into account non-sensor information, or user knowledge about the system,
in the estimation process. For instance, some state variables are known to be
nonnegative; or a control input has bounded variation; or an event is known to
happen at certain moment. All these are valuable information that can be used
for the estimation of $z$. The goal of this section is to define a measure for
the observability of $z$ using the observation data of $y$ as well as the
constraints and the control system model.

###  2.1 Definition

We assume that variables along trajectories are associated with metrics. For
instance, $y=h(t,x(t),u(t),\mu)$, as a function of $t$, has $L^{2}$ or
$L^{\infty}$ norm; $z=e(t,x(t),u(t),\mu)$ can be measured by its function
norm, or by the norm of its initial value $e(t_{0},\xi(t_{0}),u(t_{0}),\mu)$.
A metric used for $z$ is denoted by $||\cdot||_{Z}$; and $||\cdot||_{Y}$
represents the metric for $y=h(t,x,u,\mu)$. The following definition is
applicable to systems with general metrics, including $L^{p}$ and
$L^{\infty}$. Unless otherwise specified, a norm $||a||$ for $a\in\Re^{k}$ is
defined by

$(a_{1}^{2}+\cdots+a_{k}^{2})^{1/2}$

For any function $h(t)$, $t\in[t_{0},t_{f}]$, its $L^{p}$-norm is defined by

$||h||_{L^{p}}=\left(\int_{t_{0}}^{t_{f}}|h(t)|^{p}dt\right)^{1/p}$

The infinity norm is defined by

$||h||_{\infty}=\lim_{p\rightarrow\infty}\left(\int_{t_{0}}^{t_{f}}|h(t)|^{p}dt\right)^{1/p}$

which equals its essential supremum value. In this paper, a triple
$(x(t),u(t),\mu)$ represents a trajectory of (1) satisfying the differential
equations as well as the constraints. Given a positive number $\epsilon>0$ and
a nominal, or true, trajectory $(x(t),u(t),\mu)$. Define

$\begin{array}[]{rcllllllll}{\cal
E}=\left\\{(\hat{x}(t),\hat{u}(t),\hat{\mu})|\;||h(t,\hat{x}(t),\hat{u}(t),\hat{\mu})-h(t,x(t),u(t),\mu)||_{Y}\leq\epsilon\right\\}\end{array}$
(2)

The number $\epsilon$ is used as an output error bound. If
$h(\hat{x}(t),\hat{u}(t),\hat{\mu})$ stays in the $\epsilon$ neighborhood of
the nominal output $h(t,x(t),u(t),\mu)$, then we consider the trajectory
$(\hat{x}(t),\hat{u}(t),\hat{\mu})$ not distinguishable from the nominal one
using output measurement. In this case, any trajectory in $\cal E$ can be
picked by an estimation algorithm as an approximation of the true trajectory
$(x(t),u(t),\mu)$. For this reason, a trajectory in $\cal E$ is called an
estimation of $(x(t),u(t),\mu)$. Similarly,
$\hat{z}=e(t,\hat{x},\hat{u},\hat{\mu})$ is an estimation of $z=e(t,x,u,\mu)$.

###### Definition 1

Given a trajectory $(x(t),u(t),\mu)$, $t\in[t_{0},t_{1}]$. Let $\epsilon>0$ be
the output error bound. Then the number $\rho_{o}(\epsilon)$ is defined as
follows

$\begin{array}[]{lllllllll}\rho_{o}(\epsilon)=\displaystyle\max_{(\hat{x}(t),\hat{u}(t),\hat{\mu})}||e(t,\hat{x}(t),\hat{u}(t),\hat{\mu})-e(t,x(t),u(t),\mu)||_{Z}\\\
\mbox{subject to}\\\
||h(\hat{x}(t),\hat{u}(t),\hat{\mu})-h(t,x(t),u(t),\mu)||_{Y}\leq\epsilon\\\
\dot{\hat{x}}=f(t,\hat{x},\hat{u},\hat{\mu}),\\\
(\hat{x}(\cdot),\hat{u}(\cdot),\hat{\mu})\in{\cal C}\end{array}$ (3)

The number $\rho_{o}(\epsilon)$ is called the ambiguity in the estimation of
$z$ along the trajectory $(x(t),u(t),\mu)$.

Let $U$ be an open set in $(x,u,\mu)$-space and $[t_{0},t_{1}]$ be a time
interval. Then the largest value of ambiguity along all trajectories in $U$ is
called the ambiguity in the estimation of $z$ in the region $U$.

Remarks

1\. The ratio $\rho_{o}(\epsilon)/\epsilon$ measures the sensitivity of
estimation to the noise in $y$. A small sensitivity value implies strong
observability of $z$ in the presence of sensor noise.

2\. The ratio $\rho_{o}(\epsilon)/\epsilon$ is closely related to the
observability gramian. Consider a linear system

$\dot{x}=Ax,\;\;y=Cx$

Suppose $z=x$ and suppose $||\cdot||_{Y}$ is the $L^{2}$-norm. Let $P$ be the
observability gramian [1], [4], then

$\begin{array}[]{lllllllll}||y||_{Y}^{2}=x_{0}^{T}Px_{0}\end{array}$ (4)

Given $||y||_{Y}=\epsilon$, $\rho_{o}$ equals the maximum value of $||x_{0}||$
satisfying (4). In this case, $x_{0}$ is an eigenvector of $P$ associated to
the smallest eigenvalue $\lambda_{\min}$. We have

$\begin{array}[]{lllllllll}\epsilon^{2}=\lambda_{min}\rho_{o}^{2}\end{array}$
(5)

Therefore, the ratio $\rho_{o}(\epsilon)^{2}/\epsilon^{2}$ equals the
reciprocal of the smallest eigenvalue of the observability gramian. For
nonlinear systems, one can use observability gramian to approximate
$\rho_{o}(\epsilon)/\epsilon$. An advantage of this approach is that the
gramian can be computed empirically without solving the optimization problem
(3). Details on empirical computational algorithms for the gramian of
nonlinear systems can be found in [4] and [5].

3\. Given a fixed number $\rho_{o}>0$, from (4) the least sensitive initial
state defined in (5), i.e. the eigenvector of length $\rho_{o}$ associated to
$\lambda_{min}$, can be found by the following optimization

$\arg\min_{||x_{0}||_{X}=\rho_{o}}||y||_{Y}=x_{0}^{T}Px_{0}$

Extending this idea to nonlinear systems, the least observable direction in
initial states can be defined as follows: given $\rho_{o}>0$, let $\epsilon$
be the minimum value from the following problem of minimization

$\begin{array}[]{lllllllll}\epsilon=\displaystyle\min_{(\hat{x}(t),\hat{u}(t),\hat{\mu})}||h(t,\hat{x}(t),\hat{u}(t),\hat{\mu})-h(t,x(t),u(t),\mu)||_{Z}\\\
\mbox{subject to}\\\ ||\hat{x}_{0}-x_{0}||_{X}=\rho_{o}\\\
\dot{\hat{x}}=f(t,\hat{x},\hat{u},\hat{\mu}),\\\
(\hat{x}(\cdot),\hat{u}(\cdot),\hat{\mu})\in{\cal C}\end{array}$

The resulting initial $x_{0}$ represents the least observable state on the
sphere $||\hat{x}_{0}-x_{0}||_{X}=\rho_{o}$ and the ratio $\rho_{o}/\epsilon$
measures the unobservability of initial states. This definition is a reverse
process of Definition 1. However, in some cases it is easier to handle the
constraint $||\hat{x}_{0}-x_{0}||_{X}=\rho_{o}$ than the inequality of
$h(t,x,u,\mu)$ in (3).

4\. The metric for output in Definition 1 can be a vector valued function
which is bounded by a vector $\epsilon$. This flexibility is useful for
systems using different types of sensors with different accuracy.

5\. Definition 1 is independent of estimation methods. It characterizes a
fundamental attribute of the system itself, not the accuracy of a specific
estimation method. In the following, we compare Definition 1 to traditional
definitions of observability. It is shown that simple linear systems
observable in the traditional sense might be weakly observable or practically
unobservable under Definition 1; and, on the other hand, some systems not
observable under traditional definitions are practically observable with a
small ambiguity in estimation. $\diamond$

###  2.2 Computational dynamic optimization

The problem defined by (3) is a dynamic optimization. To apply Definition 1,
this problem must be solved. Obviously, an analytic solution to (3) is very
difficult to derive, if not impossible, especially in the case of nonlinear
systems. However, there exist numerical approaches that can be used to find
its approximate solution. For instance, various numerical methods are
discussed in detail in [6], [7], and [8]. Surveys on numerical methods for
solving nonlinear optimal control problems can be found in [9, 10]. The
computational algorithm used in this paper is from a family of approaches
called direct method [11, 12, 13, 14, 15]. The essential idea of this method
is to discretize the optimal control problem and then solve the resulting
finite-dimensional optimization problem. The simplicity of direct methods
makes it an ideal tool for a wide variety of applications of dynamic
optimization with constraints, including (3) in Definition 1.

More specifically, all simulations in this paper use a pseudospectral optimal
control method. In this approach, a set of nodes is selected using either the
zeros or the critical points of orthogonal polynomials, in our case the
Legendre-Guass-Lobatto nodes. Then, the problem of dynamic optimization is
discretized at the nodes to result in a nonlinear programming, which is solved
using sequential quadratic programming. Details are referred to [12, 14, 15].
In some of the following examples, dynamic optimizations are solved using the
software package DIDO [16].

A frustration in nonlinear programming is the difficulty of finding global
optimal solutions within a given domain. This is no exception in this paper.
In all examples, a variety of initial guesses are used to gain a comfortable
level of confidence that the result is not a local optimal solution. However,
for all examples of nonlinear systems in this paper, the computation cannot
guarantee global maximum value for (3). Nevertheless, in the case that a
result is not the global maximum value, it still provides a lower bound of the
ambiguity value $\rho_{o}$.

###  2.3 Examples

For the rest of this section, we illustrate Definition 1 using several
examples. In the first example, it is shown that the traditional concept of
observability is ineffective for systems with large dimensions. It justifies
the necessity of a quantitative definition of observability, such as
Definition 1.

Example. Consider the following linear system

$\begin{array}[]{rcllllllll}&&\dot{x}_{1}=x_{2}\\\ &&\dot{x}_{2}=x_{3}\\\
&&\vdots\\\
&&\dot{x}_{n}=-\displaystyle\sum_{i=1}^{n}\left(\begin{array}[]{ccc}n\\\
i-1\end{array}\right)x_{i}\\\ &&y=x_{1}\end{array}$ (6)

Under a traditional definition of observability, this system is perfectly
observable for any choice of $n$, i.e. given an output history $y=x_{1}(t)$,
it corresponds to a unique initial state $x_{0}$. However, if Definition 1 is
applied to measure the observability, it is a completely different story when
the dimension is high.

Suppose the goal is to estimate $x_{0}$. We can define $z=x(t)$. Definition 1
is applicable with arbitrary metrics. To measure the observability of the
initial state, we can use the norm of $x(0)$ as the metric for $z$, i.e.

$||z(t)||_{Z}=||x(0)||$

For this example, the output accuracy is measured by $L^{\infty}$-norm,

$||y(t)||_{Y}=\displaystyle\max_{t\in[t_{0},t_{f}]}|y(t)|$

Let us assume that the true initial state is

$x_{0}=\left[\begin{array}[]{ccccccccc}0&0&\cdots&1\end{array}\right]^{T}$

Let the output error bound be small, $\epsilon=10^{-6}$. So, we assume very
accurate observation data. The time interval is $[0,15]$. Problem (3) has the
following form

$\begin{array}[]{lllllllll}\rho_{o}(\epsilon)=\displaystyle\max_{\hat{x}}||\hat{x}(0)-x(0)||\\\
\mbox{subject to}\\\ ||\hat{x}_{1}(t)-x_{1}(t)||_{Y}\leq\epsilon\\\
\dot{\hat{x}}=f(\hat{x})\end{array}$ (7)

It is solved to compute $\rho_{o}(\epsilon)$. Table 1 lists the result for
$n=2,3,\cdots,9$.

n | 2 | 3 | 4 | 5
---|---|---|---|---
$\rho_{o}(\epsilon)$ | $4.70\times 10^{-6}$ | $2.67\times 10^{-5}$ | $1.53\times 10^{-4}$ | $8.89\times 10^{-4}$
$\epsilon$ | $10^{-6}$ | $10^{-6}$ | $10^{-6}$ | $10^{-6}$
n | 6 | 7 | 8 | 9
$\rho_{o}(\epsilon)$ | $5.20\times 10^{-3}$ | $3.01\times 10^{-2}$ | $1.75\times 10^{-1}$ | $1.02$
$\epsilon$ | $10^{-6}$ | $10^{-6}$ | $10^{-6}$ | $10^{-6}$
Table 1: Observability

From the table, when $n=2$ the ambiguity in the estimation of $x_{0}$ is as
small as $4.70\times 10^{-6}$. So the system is strongly observable.
Equivalently, if the observation data has absolute error less than $\epsilon$,
then the worst possible estimation of $x_{0}$ has an error at the scale of
$10^{-6}$. This conclusion agrees with the traditional theory of
observability. However, when the dimension is increased, the observability
ambiguity increases too; thus the system becomes less observable. At $n=8$,
the observability ambiguity is as big as $0.175$, or the worst error of
estimation is $17.5\%$ relative to the true $x_{0}$. When $n=9$, the
observability ambiguity is $1.02$. In this case, the worse relative error in
estimation is more than $100\%$! Thus, the system is practically unobservable,
although it is perfectly observable under a traditional definition. Figure 1
shows why this system is practically unobservable. The continuous curves
represent the true trajectory and its output for $n=9$; the dotted curves are
the estimation. The outputs of both trajectories agree to each other very well
(Figure on top), but the initial states (only $x_{9}$ is plotted) are
significantly different.

Figure 1: Estimation error ($n=9$)

As shown in Figure 1, while the estimation of the initial state is inaccurate,
the estimation is very close to the true value at the final time $t_{f}$. To
see the observability of the final state, let us use a different metric for
$z$,

$\begin{array}[]{lllllllll}||z(t)||_{Z}=||x(t_{f})||\end{array}$ (8)

If we consider $t=t_{f}$ as the current time moment, then this metric is used
to measure the detectability of the current system state, rather than the
observability of the initial value $x_{0}$. To compute the ambiguity under the
new metric, we solve the problem defined in (7) except that the cost function
is replaced by the metric (8). For the case of $t_{f}=10$ and
$\epsilon=10^{-6}$, the ambiguity in the estimation of $x(t_{f})$ equals
$2.7328\times 10^{-6}$. Therefore, the system is accurately detectable.

To summarize, this example shows a set of linear systems that are observable
under conventional definition. However, as the dimension is increased, the
systems become practically unobservable in the sense that an output trajectory
cannot accurately determine the state trajectory. Meanwhile, the detectability
of the system is not changed with the dimension. Definition 1 is used here to
treat both observability and detectability in the same framework,
quantitatively. $\diamondsuit$

A concept is useful only if it is verifiable for a wide spectrum of systems
and applications. An advantage of Definition 1 is that the dynamic
optimization (3) can be numerically solved for various types of applications.
In the following, we illustrate the usefulness of the ambiguity in estimation
using two examples. One is a networked cooperative control system; and the
other one is parameter identification for nonlinear systems.

Example (Partial observability of cooperative and networked systems). In this
example, it is shown that an unobservable system under traditional control
theory can be practically observable by employing user knowledge about the
system, such as an approximate upper bound of control input. Consider a
networked control system showing in Figure 2. Suppose it consists of an
unknown number of vehicles. Due to the large number of subsystems, it could be
either impossible or unnecessary to process or collect all information about
the entire system. A practical approach is to find partial observability with
local sensor information only. In this example, we assume that the cooperative
relationships in the system is unknown except that we know Vehicle 2 follows
Vehicle 1; Vehicle 3 follows both Vehicle 1 and 2 as shown by the arrows in
Figure 2. The dashed lines in the figure represent unknown cooperative
relationships. The question to be answered is the observability of Vehicle 1
if the locations of Vehicle 2 and 3 can be measured.

Figure 2: Cooperative networked system

Suppose each vehicle can be treated as a point mass with a linear dynamics

$\begin{array}[]{lllllllll}\dot{x}_{i1}=x_{i2}&\dot{y}_{i1}=y_{i2}\\\
\dot{x}_{i2}=u_{i}&\dot{y}_{i2}=v_{i}\end{array}$

Assume that the control input of vehicles 2 and 3 are defined as follows

$\begin{array}[]{lllllllll}u_{2}=a_{1}(x_{21}-x_{11}-d_{1})+a_{2}(x_{22}-x_{12})\\\
u_{3}=b_{1}(x_{31}-\displaystyle\frac{x_{11}+x_{21}}{2}-d_{2})+b_{2}(x_{32}-\displaystyle\frac{x_{12}+x_{22}}{2})\end{array}$

where $d_{i}$ is the distance of separation. The control in the $y$-direction
is the same. So, Vehicle 2 follows Vehicle 1, Vehicle 3 follows the average
position of Vehicle 1 and 2. Suppose we can measure the positions of Vehicles
2 and 3.

$\begin{array}[]{lllllllll}\mbox{output}=\left[\begin{array}[]{ccccccccc}x_{21}&y_{21}&x_{31}&y_{31}\end{array}\right]^{T}\end{array}$
(9)

The question to be answered is the observability of the location and velocity
of Vehicle 1, i.e. $x_{11}$, $y_{11}$, $x_{12}$, $y_{12}$. We would like to
emphasize that the control input of Vehicle 1 is unknown because its input is
determined by its cooperative relationships with other vehicles or agents,
which is not given. Therefore, in traditional control theory, Vehicle 1 is
unobservable.

To make Vehicle 1 practically observable with limited local measurement, we
assume that the input of Vehicle 1 has bounded variation with an upper bound
$V_{max}$. This is to say that the vehicles are not supposed to make high
frequency zigzag movement, or the control does not have chattering phenomenon.
However, discontinuity in control, such as bang-bang, is allowed. In the
following, we measure the observability of Vehicle 1 along a trajectory
defined by the parameters in Table 2.

$t_{0}$ | $t_{f}$ | $d_{1}$ | $d_{2}$ | $a_{1}$ | $a_{2}$ | $b_{1}$ | $b_{2}$ | $(x_{11}^{0},x_{12}^{0})$ | $(x_{21}^{0},x_{22}^{0})$ | $(x_{31}^{0},x_{32}^{0})$
---|---|---|---|---|---|---|---|---|---|---
$0$ | $20$ | $-2$ | $-2$ | $-1$ | $-2$ | $-3$ | $-7$ |  ($0,4$) |  ($d_{1},4$) |  ($d_{2},4$)
Table 2: Parameters of Nominal Trajectory

The control input of the nominal trajectory is

$u_{1}=\sin\displaystyle\frac{(t_{f}-t_{0})t}{\pi}$

which is unknown to the observer. To measure the ambiguity in the estimation
of vehicles, we assume the output error bound of (9) is $\epsilon=10^{-2}$.
For the unknown $u_{1}$, we assume a bounded variation of less than or equal
to $V_{max}=3.0$, which is $50\%$ higher than the true variation. The metric
for each output variable is the $L^{\infty}$-norm. The metric for the location
and velocity of Vehicle 1 is the $L^{2}$-norm. The ambiguity in the estimation
of each state variable is computed by solving a problem of dynamic
optimization. Using the estimation of $x_{11}$ as an example, we have

$\begin{array}[]{lllllllll}\rho_{o}(\epsilon)=\displaystyle\max_{(\hat{x},\hat{u}_{1})}||\hat{x}_{11}(t)-x_{11}(t)||_{L^{2}}\\\
\mbox{subject to}\\\
||\hat{x}_{21}(t)-x_{21}(t)||_{L^{\infty}}\leq\epsilon_{1}\\\
||\hat{x}_{31}(t)-x_{31}(t)||_{L^{\infty}}\leq\epsilon_{2}\\\
\dot{\hat{x}}=f(\hat{x})\\\ V(\hat{u}_{1})\leq V_{max}\end{array}$ (10)

where $V(\hat{u}_{1})$ is the total variation. In computation, this constraint
is discretized at a set of node points $t_{0}<t_{1}<\cdots<t_{N}=f_{f}$ so
that

$\displaystyle\sum_{k=1}^{N}|\hat{u}_{1}(t_{i})-\hat{u}_{1}(t_{i-1})|\leq
V_{max}$

An interesting point in the formulation (10) is that the outputs for the two
vehicles, i.e. $\hat{x}_{21}$ and $\hat{x}_{31}$, have different error bounds,
$\epsilon_{1}$ and $\epsilon_{2}$. The metric for the outputs is a vector
valued function. This flexibility of using different $\epsilon$ value for
multiple outputs is advantageous for systems with multiple sensors of
different qualities that measure various states with different accuracy.

The computed result is shown in Table 3. The small relative ambiguity value
shows that the location and velocity of Vehicle 1 are practically observable
given the measurement of the positions of Vehicle 2 and 3, without using any
information about the rest of the networked system and without knowing the
input of Vehicle 1. The worst estimation of $x_{11}$ and $x_{12}$ is shown in
Figure 3, which has good accuracy. $\diamondsuit$

$V_{max}$ | $\epsilon$ | $\rho_{x_{11}}$ | $\rho_{x_{11}}/||x_{11}||_{L_{2}}$ | $\rho_{x_{12}}$ | $\rho_{x_{12}}/||x_{12}||_{L_{2}}$
---|---|---|---|---|---
$3$ | $10^{-2}$ | $1.2257$ | $2.8\times 10^{-3}$ | $0.5901$ | $1.16\times 10^{-2}$
Table 3: Observability of Vehicle 1 Figure 3: The worst estimation of the
position and velocity of Vehicle 1

In the following example, the concept of ambiguity in estimation is applied to
the Laub-Loomis model [17] with unknown parameters, a nonlinear system of
oscillating biochemical network.

Example (Parameter identification) In the study of biochemical networks, it
was proposed that interacting proteins could account for the spontaneous
oscillations in adenylyl cyclase activity that was observed in homogeneous
populations of dictyostelium cells. While a set of terminologies such as
$3^{\prime},5^{\prime}$-cycle monophosphate (cAMP) and adenylate cyclase (ACA)
are involved in the problem, we focus on the state space in which a set of
seven nonlinear differential equations are used as the model [17].

$\begin{array}[]{rcllllllll}\dot{x}_{1}&=&k_{1}x_{7}-k_{2}x_{1}x_{2}\\\
\dot{x}_{2}&=&k_{3}x_{5}-k_{4}x_{2}\\\
\dot{x}_{3}&=&k_{5}x_{7}-k_{6}x_{2}x_{3}\\\
\dot{x}_{4}&=&k_{7}-k_{8}x_{3}x_{4}\\\
\dot{x}_{5}&=&k_{9}x_{1}-k_{10}x_{4}x_{5}\\\
\dot{x}_{6}&=&k_{11}x_{1}-k_{12}x_{6}\\\
\dot{x}_{7}&=&k_{13}x_{6}-k_{14}x_{7}\\\ \end{array}$ (11)

In a robustness study [18], it was shown that a small variation in the model
parameters can effectively destroy the required oscillatory dynamics. As a
related question, it becomes interesting to investigate the possibility of
estimating the parameters in the system, $k_{1}$, $k_{2}$, $\cdots$, $k_{7}$.
To exemplify the idea, we assume that $x_{1}$, the value of CAC, is
measurable, i.e.

$y=x_{1}$

We also assume that the initial states in experimentation is known. Suppose
the unknown parameters are $k_{1}$, $k_{6}$, and $k_{10}$; and suppose the
other parameters are known. The goal is to use the measured data of $y$ to
estimate the unknown parameters , i.e.

$z=\left[\begin{array}[]{ccccccccc}k_{1}&k_{6}&k_{10}\end{array}\right]$

Using Definition 1 we can quantitatively determine the observability of the
unknown parameters. Along a nominal trajectory $(x^{\ast}(t),k^{\ast})$, the
ambiguity can be computed by solving the following special form of (3).

$\begin{array}[]{lllllllll}\rho_{o}^{2}=\displaystyle{\max_{(x,k_{1},k_{6},k_{10})}}(k_{1}-k^{\ast}_{1})^{2}+(k_{6}-k^{\ast}_{6})^{2}+(k_{10}-k^{\ast}_{10})^{2}\\\
\mbox{subject to}\\\
||x_{1}(t)-x^{\ast}_{1}(t)||_{L^{2}}^{2}\leq\epsilon^{2}\\\
\dot{x}=f(t,x,k_{1},k_{6},k_{10}),\mbox{ other parameters equal nominal
value}\\\ x(t_{0})=x^{\ast}(t_{0})\\\ \end{array}$

In the simulation, a nominal trajectory is generated using the following
parameter value and initial condition

$\begin{array}[]{lllllllll}k_{1}=2.0,\;\;k_{2}=0.9,\;\;k_{3}=2.5,\;\;k_{4}=1.5,\;\;k_{5}=0.6,\;\;k_{6}=0.8,\;\;k_{7}=1.0,\\\
k_{8}=1.3,\;\;k_{9}=0.3,\;\;k_{10}=0.8,\;\;k_{11}=0.7,\;\;k_{12}=4.9,\;\;k_{13}=23.0,\;\;k_{14}=4.5,\\\
x(0)=\left[\begin{array}[]{ccccccccc}1.9675&1.2822&0.6594&1.1967&0.6712&0.2711&1.3428\end{array}\right]\end{array}$

The output error bound is being set at $\epsilon=10^{-2}$. The computation
reveals that the ambiguity in the estimation of
$z=\left[\begin{array}[]{ccccccccc}k_{1}&k_{6}&k_{10}\end{array}\right]$ is

$\rho_{o}=2.38\times 10^{-2}$

Given the nominal value of the parameters, the relative ambiguity in
estimation is about $1\%$. So, the parameters are strongly observable. In
fact, the worst estimation of the parameters is

$k_{1}=2.0150,\;\;k_{6}=0.8082,\;\;k_{10}=0.7836.$

The trajectory generated by the worst parameter estimation is shown in Figure
4. $\diamond$

Figure 4: The trajectory of worst estimation: curve - true trajectory; star -
estimation)

##  3 Input-to-output gain

$L^{p}$-gain is a tool of analysis widely used by control engineers to
quantitatively measure the sensitivity and robustness of systems. Consider

$\begin{array}[]{rcllllllll}\dot{x}=f(t,x,w,\mu)\\\ z=e(t,x,w,\mu)\end{array}$
(12)

where $x\in\Re^{n_{x}}$ is the state variable, $w\in\Re^{n_{w}}$ is the input
that represents the disturbance, $\mu\in\Re^{n_{\mu}}$ is the system
uncertainty or a parameter, $z\in\Re^{n_{z}}$ is the performance. In the
following, the $L^{p}$-norm of a vector valued function is denoted by
$||\cdot||_{L^{p}}$, for instance

$||z(t)||_{L^{p}}=\displaystyle\left(\int_{t_{0}}^{t_{1}}\sum_{i=1}^{n_{z}}|z_{i}(t)|^{p}dt\right)^{1/p}$

Given a fixed time interval $[t_{0},t_{1}]$ and $\sigma>0$. Suppose the input
$w$ is a function in $L^{p}$ space for some $1\leq p\leq\infty$ such that
$w(t)$ is bounded by $\sigma$, i.e.

$||w(t)||_{L^{p}}\leq\sigma$

Let $x^{\ast}(t)$ be a nominal trajectory with $x^{\ast}(t_{0})=x_{0}$,
$w^{\ast}(t)=0$. Fix the initial value $x(t_{0})=x_{0}$. Suppose the system
uncertainty is bounded, $\mu_{min}\leq\mu\leq\mu_{max}$. Then the $L^{p}$-gain
from $w$ to $z$ along $x^{\ast}(t)$ is defined as follows

$\begin{array}[]{lllllllll}\gamma(\sigma)=\displaystyle\max_{\begin{array}[]{c}||w||_{L^{p}}\leq\sigma,\\\
\mu_{min}\leq\mu\leq\mu_{max}\end{array}}\displaystyle\frac{||e(t,x,w,\mu)-e(t,x^{\ast},0,\mu)||_{L^{p}}}{\sigma}\end{array}$
(13)

Remark 6. Without parameter, the maximum value of
$||e(t,x,w)-e(t,x^{\ast},0)||_{L^{p}}$ is the ambiguity in the estimation of
$z=e(t,x,w)$. More specifically, consider the ambiguity in the estimation of
$z$ under the observation of $w$ with an error bound $\sigma$. Then the
$L^{p}$-gain $\gamma(\sigma)$ gain equals the ratio of the ambiguity and
$\sigma$.

###  3.1 Computation and example

The input-to-output gain can be computed by solving the problem (13). In
Section 2, the output function $y$ is smooth. It can be numerically
approximated in a finite dimensional space, such as interpolation using a
finite number of nodes. Similarly, one has to work on a finite dimensional
space of $w$ to carry out the computation. So, for the purpose of computation,
we only discuss the $L^{p}$-gain in a finite dimensional space of $w$, denoted
by $\cal U$, rather than the infinite dimensional space of arbitrary
integrable functions. The space ${\cal U}$ can be defined by the frequency
bandwidth, or the order of polynomials, or some other spaces used for the
approximation of the input. Then, (13) is reformulated as follows

$\begin{array}[]{lllllllll}\gamma_{\cal
U}(\sigma)=\displaystyle\max_{\begin{array}[]{c}w\in{\cal
U},||w||_{L^{p}}\leq\sigma\\\
\mu_{min}\leq\mu\leq\mu_{max}\end{array}}\displaystyle\frac{||e(t,x,w,\mu)-e(t,x^{\ast},0,\mu)||_{L^{p}}}{\sigma}\end{array}$
(14)

More specifically, given a positive number $\sigma>0$ define

$J(x(\cdot),w(\cdot),\mu)=||e(t,x(\cdot),w(\cdot),\mu)-e(t,x^{\ast}(\cdot),0,\mu)||_{L^{p}}$

Then the following dynamic optimization determines the $L^{p}$-gain over the
space ${\cal U}$.

Dynamic optimization for $L^{p}$-gain

$\begin{array}[]{lllllllll}\rho=\displaystyle{\max_{(x,w,\mu)}}J\\\
\mbox{subject to}\\\ w(t)\in{\cal U},\;||w(t)||_{L^{p}}\leq\sigma\\\
\dot{x}=f(t,x,w,\mu),\\\ x(t_{0})=x_{0}\\\
\mu_{min}\leq\mu\leq\mu_{max}\end{array}$ (15)

The $L^{p}$-gain from $w\in{\cal U}$ to $z$ is $\gamma_{\cal
U}(\sigma)=\displaystyle\frac{\rho}{\sigma}$. $\diamond$

Example ($L^{p}$-gain in the presence of system uncertainty). Atomic force
microscope (AFM) invented two decades ago is used to probe surfaces at the
atomic level with good accuracy. This type of equipment is also used as nano-
manipulation tools to handle particles at nano-scale [19, 20, 21]. Illustrated
in Figure 5, the system consists of a microcantilever with a sharp tip at one
end. The vibration of the cantilever is measured by an optical sensor. The
topographic images of surfaces can be taken by measuring the cantilever’s
dynamic behavior which is determined by the interacting force of the tip with
the sample.

Figure 5: Atomic force microscope

The dynamics of the vibrating tip can be modeled as a second order system [21]

$\begin{array}[]{lllllllll}\dot{x}_{1}=x_{2}\\\
\dot{x}_{2}=-\omega^{2}x_{1}-2\xi\omega x_{2}+h(x_{1},\delta)+u(t)+w(t)\\\
z=x_{1}\end{array}$ (16)

where $x_{1}$ is the position of the cantilever tip at the scale of
nanometers, $x_{2}$ is its velocity, $\omega$ is the natural frequency of the
cantilever, and $\xi$ is the damping coefficient. In this system, $u(t)$ is
the control input, $w(t)$ is the actuator disturbance which is unknown. The
function $h(x_{1},\delta)$ is the tip-sample interaction force in which
$\delta$ is the separation between the equilibrium of $x_{1}$ and the sample
surface. It is a system uncertainty. We adopt the following model for
$h(x_{1},\delta)$ [20].

$h(x_{1},\delta)=-\displaystyle\frac{\alpha_{1}}{(\delta+x_{1})^{2}}+\displaystyle\frac{\alpha_{2}}{(\delta+x_{1})^{8}}$

At nano-scale, system uncertainly and performance robustness are critical
issues in control design because a seemingly small noise or uncertainty may
have significant impact on the performance. In the following, we assume that
the value of $\delta$ and $w(t)$ are unknown. The goal is to compute the
$L^{2}$-gain from the actuator disturbance $w$ to the performance $z=x_{1}$
for $\delta$ in the entire interval $[\delta_{min},\delta_{max}]$.

The simulations are based on the following set of parameter value

$\begin{array}[]{lllllllll}\omega=1.0,&\xi=0.02,&\alpha_{1}=0.1481,&\alpha_{2}=3.6\times
10^{-6}\end{array}$

The nominal control input is $u(t)=1$ and the time interval is $[0,7]$, which
is long enough to cover one period of oscillation. The bounds of $\delta$ are

$\begin{array}[]{lllllllll}\delta_{min}=0.8,&\delta_{max}=1.2\end{array}$

We use $L^{2}$-norm as the metric for the actuator disturbance force $w$. Let
$w$ be an arbitrary function in a two-frequency space ${\cal W}_{k_{1},k_{2}}$
defined as follows

$\begin{array}[]{lllllllll}w=\displaystyle\sum_{i=1}^{2}\left(A_{i}\cos(\displaystyle\frac{2\pi
k_{i}}{t_{f}-t_{0}}t)+B_{i}\sin(\displaystyle\frac{2\pi
k_{i}}{t_{f}-t_{0}}t)\right)\end{array}$

Let $\sigma=0.03$. Then the $L^{2}$-gain is computed by solving the dynamic
optimization defined in (15). More specifically,

$\begin{array}[]{lllllllll}\rho=\displaystyle{\max_{(x,u,\delta)}}||x_{1}(t)-x_{1}^{\ast}(t)||_{L^{2}}\\\
\mbox{subject to}\\\ ||w(t)||_{L^{2}}\leq\sigma,\;w(t)\in{\cal
W}_{k_{1},k_{2}}\\\ \dot{x}=f(t,x,u,\delta),\\\ \hat{x}(t_{0})=x_{0}\\\
\delta\in[\delta_{min},\delta_{max}]\end{array}$

The $L^{2}$-gain equals $\displaystyle\frac{\rho}{\sigma}$. It is computed for
spaces with various frequencies, ${\cal W}_{0,1}$, ${\cal W}_{2,3}$, ${\cal
W}_{4,5}$, ${\cal W}_{6,7}$, and ${\cal W}_{8,9}$. The result is shown in
Figure 6. The $L^{2}$-gain for frequencies $0$ and $1$ is $2.5707$. When the
frequencies are increased, the gain decreases. $\diamond$

Figure 6: $L^{2}$-gain

###  3.2 An alternative algorithm for $L^{2}$-gain

Solving (15) becomes increasingly difficult for high frequencies. The reason
is that, for the computational purpose, the problem of dynamic optimization is
always discretized at finite number of nodes in time. For higher frequencies,
the number of nodes must be increased. As a result, the dimension of
optimization variables is increased as well. Developing efficient methods of
computation for inputs of high frequencies requires further research. However,
in the case of $L^{2}$-gain, there exists an alternate approach without
solving dynamic optimization. Inspired by Remark 6 and the observability
gramian in [4], the $L^{2}$-gain can be approximated using the following
matrix.

Suppose the space of input, ${\cal U}$, is finite dimensional with a basis,
$w_{1}$, $w_{2}$, $\cdots$, $w_{m}$. Let $\sigma>0$ be a constant number. For
the input $\pm\sigma w_{i}$, the trajectory of

$\begin{array}[]{lllllllll}\dot{x}=f(t,x,\pm\sigma w_{i})\\\
x(t_{0})=x_{0}\end{array}$

is denoted by $x^{i\pm}(t)$. Define

$\begin{array}[]{lllllllll}2\Delta z^{i}=e(t,x^{i+}(t),\sigma
w_{i}(t))-e(t,x^{i-}(t),-\sigma w_{i}(t))\end{array}$ (17)

Now, define

$\begin{array}[]{lllllllll}G^{w}_{ij}=<w_{i},w_{j}>=\displaystyle\frac{1}{t_{f}-t_{0}}\displaystyle\int_{t_{0}}^{t_{f}}w_{i}(t)^{T}w_{j}(t)dt\\\
G^{z}_{ij}=<\Delta z^{i},\Delta
z^{j}>=\displaystyle\frac{1}{t_{f}-t_{0}}\displaystyle\int_{t_{0}}^{t_{f}}\Delta
z^{i}(t)^{T}\Delta z^{j}(t)dt\end{array}$ (18)

Denote the matrices $G^{w}=(G^{w}_{ij})_{i,j=1}^{n}$ and
$G^{z}=(G^{z}_{ij})_{i,j=1}^{n}$. Given any

$w=\displaystyle\sum_{i=1}^{m}a_{i}w_{i}$

satisfying

$<w,w>=\sigma^{2}$

Then $\Delta z$ generated by $\pm w$ is approximately

$2\Delta z=e(t,x^{+}(t),w(t))-e(t,x^{-}(t),-w(t))\approx
2\displaystyle\sum_{i=1}^{m}a_{i}\Delta z^{i}$

Therefore,

$\begin{array}[]{rcllllllll}||w||_{L^{2}}^{2}&=&\left[\begin{array}[]{ccccccccc}a_{1}\cdots
a_{m}\end{array}\right]G^{w}\left[\begin{array}[]{ccccccccc}a_{1}&\cdots&a_{m}\end{array}\right]^{T}\\\
||\Delta z||_{L^{2}}^{2}&\approx&\left[\begin{array}[]{ccccccccc}a_{1}\cdots
a_{m}\end{array}\right]G^{z}\left[\begin{array}[]{ccccccccc}a_{1}&\cdots&a_{m}\end{array}\right]^{T}\end{array}$

Therefore, the $L^{p}$-gain square is approximately the solution of the
following optimization

$\begin{array}[]{lllllllll}\displaystyle\frac{1}{\sigma^{2}}\max_{a}a^{T}G^{z}a\\\
\mbox{ subject to}\\\ a^{T}G^{w}a=\sigma^{2}\end{array}$

where
$a=\left[\begin{array}[]{ccccccccc}a_{1}&a_{2}&\cdots&a_{m}\end{array}\right]^{T}$.
A necessary condition for the optimal solution is

$G^{z}a=\lambda G^{w}a$

for some scalar $\lambda$. At this point,

$\begin{array}[]{rcllllllll}a^{T}G^{z}a&=&\lambda a^{T}G^{w}a\\\
&=&\lambda\sigma^{2}\end{array}$

Therefore,

$\begin{array}[]{rcllllllll}\gamma_{\cal
U}(\sigma)^{2}&=&\displaystyle\max_{||w||_{L^{2}}=\sigma^{2}}\displaystyle\frac{||\Delta
z||_{L^{2}}^{2}}{||w||_{L^{2}}^{2}}\\\
&\approx&\displaystyle\frac{a^{T}G^{z}a}{\sigma^{2}}\\\ &=&\lambda\end{array}$

On the other hand, $\lambda$ is an eigenvalue of $(G^{w})^{-1}G^{z}$. So, the
$L^{2}$-gain is approximately the square root of the largest eigenvalue.

To summarize, given a system

$\begin{array}[]{lllllllll}\dot{x}=f(t,x,w)\\\ x(t_{0})=x_{0}\end{array}$

and a space of input functions ${\cal U}$ with basis $w_{1}$, $w_{2}$,
$\cdots$, $w_{m}$. Given $\sigma>0$, compute $\Delta z^{i}$ in (17). Compute
the matrices $G^{w}$ and $G^{z}$ in (18). Then the $L^{2}$-gain is
approximately $\sqrt{\lambda_{max}}$, where $\lambda_{max}$ is the largest
eigenvalue of $(G^{w})^{-1}G^{z}$.

Example. Consider the model of AFM defined in (16). Assume that the value of
$\delta=1.0$ is known. The other parameters are the same as in the previous
example. The approximate $L^{2}$-gain is computed using the matrix approach.
The result is shown in Figure 7. This method is straightforward in computation
because no optimization is required. However, the approximation does not take
into full account the nonlinear dynamics. Comparing to the gain using (15),
the errors of the gain computed using covariance matrix is around $18\sim
20\%$ in ${\cal W}_{0,1}$, ${\cal W}_{2,3}$, and ${\cal W}_{4,5}$.

Figure 7: Approximate $L^{2}$-gain

##  4 Reachability

Dynamic optimization can be applied to quantitatively measure reachability.
Consider a control system

$\begin{array}[]{rcllllllll}\dot{x}=f(t,x,u)\end{array}$ (19)

where $x\in\Re^{n_{x}}$ and $u\in\Re^{n_{u}}$. Suppose $||x||_{X}$ is a norm
in $\Re^{n_{x}}$. We suppose that the state and control are subject to
constraint

$(x(\cdot),u(\cdot))\in{\cal C}$

###### Definition 2

Given $x_{0}$ and $x_{1}$ in $\Re^{n_{x}}$. Define

$\begin{array}[]{lllllllll}\rho_{c}(x_{0},x_{1})^{2}=\displaystyle\min_{(x,u)}||x(t_{1})-x_{1}||^{2}_{X}\\\
\mbox{subject to}\\\ \dot{x}=f(x,u)\\\ x(t_{0})=x_{0}\\\
(x(\cdot),u(\cdot))\in\cal C\end{array}$ (20)

The number $\rho_{c}(x_{0},x_{1})$ is called the ambiguity in control.

Let $D_{0},D_{1}\subset\Re^{n_{x}}$ be subsets in state space. The ambiguity
in control over the region $\bar{D}_{0}\times\bar{D}_{1}$ is defined by the
following max-min problem.

$\begin{array}[]{lll}\rho_{c}=\displaystyle{\max_{(x_{0},x_{1})\in\bar{D}_{0}\times\bar{D}_{1}}}\rho_{c}(x_{0},x_{1})\end{array}$

In this definition, $t_{1}$ is either fixed or free in a time interval
$[t_{0},T]$. In the following discussion, we assume $t_{1}$ is fixed. If
$\rho_{c}(x_{0},z_{1})$ is nonzero, then the state cannot reach $x_{1}$ from
$x_{0}$ by using admissible controls. The maximum value of
$\rho_{c}(x_{0},z_{1})$ over $\bar{D}_{0}\times\bar{D}_{1}$ represents the
worst scenario of reachability. In some applications, the relative ambiguity

$\displaystyle\frac{\rho_{c}(x_{0},x_{1})}{||x_{1}||_{X}}$

is used to measure the reachability.

The definition of $\rho_{c}$ is consistent with the classic definition of
controllability for linear time-invariant systems. To be more specific,
consider a linear system

$\begin{array}[]{rcllllllll}\dot{x}=Ax+Bu\end{array}$ (21)

let $u(\cdot)\in{\cal C}$ be the space of continuous functions from
$[t_{0},t_{1}]$ to $\Re^{m}$; and let $D_{0}=D_{1}=\Re^{n}$. If $(A,B)$ is
controllable, i.e.

$rank(\left[\begin{array}[]{ccccccccc}B&AB&A^{2}B&\cdots&A^{n-1}B\end{array}\right]=n$

then for any $x_{0}$ and $x_{1}$, there always exists a control input so that
$x(t)$ with $x(0)=x_{0}$ reaches $x_{1}$ at $t=t_{1}$. Therefore,
$\rho_{c}(x_{0},x_{1})$ is always zero for arbitrary $(x_{0},x_{1})$; and
$\rho_{c}=0$. On the other hand, if $(A,B)$ is uncontrollable, then under a
change of coordinates an uncontrollable subsystem can be decoupled from the
controllable part of the system. In the uncontrollable subsystem, the states
cannot be driven to close to each other by control inputs. Therefore,
$\rho_{c}(x_{0},x_{1})$ is unbounded for arbitrary states in the
uncontrollable subspace. This implies $\rho_{c}=\infty$. To summarize,

$\rho_{c}=\left\\{\begin{array}[]{lll}0&\mbox{if }(A,B)\mbox{ is controllable
}\\\ \infty&\mbox{if }(A,B)\mbox{ is uncontrollable }\end{array}\right.$

Control has a cost. For weakly reachable systems, it takes relatively large
control energy to reach a terminal state. The cost in reachability can be
measured by the following quantity. Denote $||u(\cdot)||_{\cal U}$ and
$||x||_{X}$ the metrics of the control input and the state, respective.

###### Definition 3

Given initial and final states, $x_{0}$ and $x_{1}$, define

$\begin{array}[]{lllllllll}W(x_{0},x_{1})=\displaystyle\min_{(x,u)}\displaystyle\lim_{\psi\rightarrow\infty}\left(||u(t)-u^{\ast}(t)||_{\cal
U}+\psi||x(t_{1})-x_{1}||_{X}\right)\\\ \mbox{subject to}\\\ \dot{x}=f(x,u)\\\
x(t_{0})=x_{0},\\\ (x(\cdot),u(\cdot))\in\cal C\end{array}$ (22)

This definition has the following property. If $x_{1}$ can be reached from
$x_{0}$, then $W(x_{0},x_{1})$ equals the minimum control

$\begin{array}[]{lllllllll}W(x_{0},x_{1})=\displaystyle\min_{(x,u)}||u(t)-u^{\ast}(t)||_{\cal
U}\\\ \mbox{subject to}\\\ \dot{x}=f(x,u)\\\
x(t_{0})=x_{0},\;\;x(t_{1})=x_{1}\\\ (x(\cdot),u(\cdot))\in\cal C\end{array}$

If $x(t)$ cannot reach $x_{1}$ using admissible control, then
$W(x_{0},x_{1})=\infty$. A large value of $W(x_{0},x_{1})$ implies higher
control cost, thus weak reachability.

Remark 7. Suppose the system is linear. Suppose $D_{0}=\\{0\\}$ and
$D_{1}=\\{x|\;\;||x||<\epsilon\\}$ for some small $\epsilon>0$. Let $W$ be the
maximum cost in reachability under $L^{2}$-norm,

$W=\displaystyle\max_{x_{1}\in\bar{D}_{1}}W(0,x_{1})$

Then $(W/\epsilon)^{2}$ equals the reciprocal of the smallest eigenvalue of
the controllability gramian.

To justify Remark 7, consider

$\dot{x}=Ax+Bu$

If $(A,B)$ is uncontrollable, we know $W=\infty$. We also know that the
smallest eigenvalue of the controllability gramian is zero. Therefore the
claim holds true. Now suppose $(A,B)$ is controllable. Define

$||u||_{\cal U}^{2}=\int_{t_{0}}^{t_{1}}||u(t)||^{2}dt$

The control cost to reach $x_{1}$ is defined by

$\begin{array}[]{lllllllll}W(0,x_{1})=\displaystyle\min_{(x,u)}||u(t)||_{\cal
U}\\\ \mbox{subject to}\\\ \dot{x}=Ax+Bu\\\
x(t_{0})=0,x(t_{1})=x_{1}\end{array}$

Let $P$ be the controllability gramian, then it is known [22] that the optimal
cost satisfies

$\begin{array}[]{rcllllllll}&&\displaystyle\int_{t_{0}}^{t_{1}}||u(t)||^{2}dt\\\
&=&x_{1}^{T}P^{-1}x_{1}\end{array}$

If $\sigma_{min}$ is the smallest eigenvalue of $P$, then

$\begin{array}[]{rcllllllll}(W/\epsilon)^{2}&&=&&\displaystyle\frac{1}{\epsilon^{2}}\max_{||x_{1}||\leq\epsilon}x_{1}^{T}P^{-1}x_{1}\\\
&&=&&\displaystyle\frac{1}{\sigma_{\min}}\end{array}$

$\diamond$

###  4.1 Example

In the following, we compute the ambiguity in the control of a heat equation
with Neumann boundary control.

$\begin{array}[]{rcllllllll}&&\displaystyle\frac{\partial w(r,t)}{\partial
t}-\kappa\displaystyle\frac{\partial^{2}w(r,t)}{\partial r^{2}}=0\\\
&&w(r,0)=0,\;\;\;0\leq r\leq 2\pi\\\ &&w(0,t)=0,\;\;\;0\leq t\leq t_{f}\\\
&&w_{rt}(r,t)|_{r=2\pi}=u(t)\end{array}$ (23)

where $w(r,t)\in\Re$ is the state of the system, $r\in\Re$ is the space
variable, and $t$ is time. The control input is $u$. For a thermal problem,
$u$ represents the rate of the heat flux $w_{r}$. The initial state is assumed
to be zero.

For the purpose of computation, we discretize the problem at equally spaced
nodes,

$0=r_{0}<r_{1}<r_{2}<\cdots<r_{N}=2\pi$

Define

$x_{1}(t)=w(r_{1},t),\;x_{2}(t)=w(r_{2},t),\;\cdots,x_{N}(t)=w(r_{N},t)$

Using central difference in space, (LABEL:pdemodel) is approximated by the
following control system defined by ODEs.

$\begin{array}[]{rcllllllll}\dot{x}_{1}&=&\kappa\displaystyle\frac{x_{2}-2x_{1}}{\Delta
r^{2}}\\\ \dot{x}_{2}&=&\kappa\displaystyle\frac{x_{1}+x_{3}-2x_{2}}{\Delta
r^{2}},\\\ &\vdots&\\\
\dot{x}_{i}&=&\kappa\displaystyle\frac{x_{i-1}+x_{i+1}-2x_{i}}{\Delta
r^{2}},\\\ &\vdots\\\
\dot{x}_{N-1}&=&\kappa\displaystyle\frac{x_{N-2}+x_{N}-2x_{N-1}}{\Delta
r^{2}}\\\ \dot{x}_{N}&=&v\\\ v&=&\dot{x}_{N-1}+\Delta ru\end{array}$ (24)

We understand that more sophisticated algorithms of solving the heat equation
exit. We adopt this central difference method for simplicity in the
illustration of control ambiguity. System (24) is linear and controllable. So,
it is theoretically a reachable system. However, in reality the maximum
temperature cannot exceed safety margin. Under such constraint, a controllable
linear system may not be reachable due to overshot. Let $t_{f}=150$ and
$\kappa=0.14$. Suppose the target states are the following arches. A few of
them is shown in Figure 8.

$w_{f}(r)=w(r,t_{f})=A\sin(r/2),\;\;\;0\leq r\leq 2\pi$

Figure 8: The target state $w_{f}(r)$

The goal is to compute the ambiguity in control from $w(r,0)=0$ to $w_{f}(r)$
with the magnitude of $0\leq A\leq 1.2$ subject to the constraint

$w(r,t)\leq 2$

The norm in the finite dimensional state space is defined by

$||x||^{2}=\Delta r\sum_{i=1}^{N}x_{i}^{2}$

which approximates the $L^{2}$-norm in $C[0,2\pi]$. To compute the control
ambiguity in the range of $0\leq A\leq 1.2$, we consider the following nodes
in the magnitude of the target arch

$A=0,\;0.2,\;0.4,\;0.6,\;0.8,\;1.0,\;1.2$

The corresponding dynamic optimization problem for the control ambiguity is
defined by

$\begin{array}[]{lllllllll}\rho_{c}(0,x_{f})^{2}=\displaystyle\min_{(x,u)}\Delta
r\sum_{i=1}^{N}(x_{i}-w_{f}(r_{i}))^{2}\\\ \mbox{subject to}\\\
\dot{x}=f(x,v)\\\ x_{i}(t)\leq 2,\;\;\;i=1,2,\cdots,N\\\ x(0)=0\end{array}$
(25)

where $f(x,v)$ is the dynamics defined in (24) with the input $v$, and

$x_{f}=\left[\begin{array}[]{ccccccccc}w_{f}(r_{1})&w_{f}(r_{2})&\cdots&w_{f}(r_{N})\end{array}\right]$

In the simulation, $N$ is selected to be $N=31$. Problem (25) is solved using
Pseudospectral method at Legendre-Gauss-Lebato (LGL) nodes [12, 14, 15]. We
use $15$ LGL nodes in this example. Through computation, it is found that the
system becomes increasingly unreachable due to the constraint when the value
of $A$ is bigger than $0.4$. The relative ambiguity in control is shown in
Figure 9. When the magnitude of the target state is $1.2$, the relative
ambiguity shows that the closest state that the system can reach has almost a
$40\%$ relative error. Therefore, the system is practically unreachable if the
state is required to be bounded.

Figure 9: Magnitude of target state vs. relative control ambiguity

##  5 Conclusion

It is shown by numerous examples and definitions that computational dynamic
optimization is a promising tool of quantitatively analyzing control system
properties. Using computational approaches, the concepts studied in this
paper, including the ambiguity in estimation and control, input-to-output
gain, and the cost in reachability, are applicable to a wide spectrum of
applications. In addition, these concepts are defined and applied in a way so
that one can take advantage of user knowledge or take into account system
constraints. As a result, the properties of control systems are not only
verified, but also measured quantitatively. While these concepts can be
applied to a wide spectrum of problems, some specific applications exemplified
in this study include: strongly observable (detectable) or weakly observable
(detectable) systems; improving observability by employing user knowledge;
partial observability of networked complex systems; $L^{2}$-gain of nonlinear
control systems; reachability in the presence of state constraints; and
boundary control of partial differential equations.

Similar to many nonlinear optimization problems, a main drawback of the
approach is that a global optimization is, in general, not guaranteed for
nonlinear systems. In addition, the problem of computational accuracy also
poses many questions remain to be answered.

## References

  * [1] T. Kailath, Linear Systems, Prentice-Hall, Inc., Englewood Cliffs, N.J., 1980.
  * [2] A. Isidori, Nonlinear Control Systems, Springer-Verlag, London, 1995.
  * [3] K. Zhou, J. Doyle, and K. Glover , Robust and Optimal Control, Prentice Hall, 1995.
  * [4] A. J. Krener and K. Ide, Measures of Unobservability, preprint, 2009.
  * [5] S. Lall, J. E. Marsden and S. Glavski, A subspace approach to balanced truncation for model reduction of nonlinear control systems, Int. J. Robust Nonlinear Control, vol. 12, 2002, pp. 519-535.
  * [6] A. E. Bryson and Y. C. Ho, Applied Optimal Control. Hemisphere, New York, 1975.
  * [7] A. E. Bryson, Dynamic Optimization, Addison-Wesley Longman, Inc., 1999.
  * [8] E. Polak, Optimization: Algorithms and Consistent Approximations, Springer-Verlag, Heidelberg, 1997.
  * [9] J. T. Betts, “Survey of Numerical Methods for Trajectory Optimization,” Journal of Guidance, Control, and Dynamics, Vol. 21, No. 2, 1998, pp. 193-207.
  * [10] E. Polak, A historical survey of computations methods in optimal control, SIAM Review, Vol. 15, pp. 553-548, 1973.
  * [11] J. T. Betts, Practical Methods for Optimal Control Using Nonlinear Programming, SIAM, Philadelphia, PA, 2001.
  * [12] G. Elnagar, M. A. Kazemi and M. Razzaghi, The pseudospectral Legendre method for discretizing optimal control problems, IEEE Trans. Automat. Contr. Vol. 40, pp. 1793-1796, 1995.
  * [13] W. W. Hager, Runge-Kutta methods in optimal control and the transformed adjoint system, Numerische Mathematik, Vol. 87, pp. 247-282, 2000.
  * [14] F. Fahroo, I. M. Ross (1998) Costate Estimation by a Legendre Pseudospectral Method. Proceedings of the AIAA Guidance, Navigation and Control Conference, 10-12 August 1998, Boston, MA.
  * [15] W. Kang, Q. Gong, and I. M. Ross, On the Convergence of Nonlinear Optimal Control using Pseudospectral Methods for Feedback Linearizable Systems, International Journal of Robust and Nonlinear Control, Vol. 17, 1251-1277, 2007.
  * [16] DIDO, Elissar LLC, http://www.elissar.biz/.
  * [17] M. T. Laub, and W. F. Loomis, A molecular network that produces spontaneous oscillations in excitable cells of Dictyostelium, Mol. Biol. Cell, Vol. 9, 3521-3532, 1998.
  * [18] J. Kim, D. G. Bates, I. Postlethwaite, L. Ma, and P. A. Iglesias, Robustness analysis of biochemical network models, IEE Proc. Syst. Biol. Vol. 153, No. 3, May, 2006.
  * [19] M. Ashhab, M.V. Salapaka, M. Dahleh, and I. Mezić, Dynamical analysis and control of microcantilevers, Automatica, 1663-1670, 1999.
  * [20] M. Basso, L. Giarre, M. Dahleh, and I. Mezle, Complex dynamics in harmonically excited Lennard-Jones oscillator : Microcantilever-sample interaction in scanning probe microscopes, ASME Journal of Dynamics Systems, Measurment and Control, Vol. 122, pp. 240-245, 2000.
  * [21] A. Delnavaz, N. Jalili, and Hassan Zohoor, Vibration control of AFM tip for nano-manipulation using combined sliding mode techniques, IEEE Proc 7th International Conference on Nanotechnology, Hong Kong, August 2 - 5, 2007.
  * [22] J. Zabczyk, Mathematical Control Theory: An Introduction, Birkhaüser, Boston, 1992.


M.~R. Hestenes and E.~Stiefel.
\newblock Methods of conjugate gradients for solving.
\newblock {\em Journal of research of the National Bureau of Standards},
49(6):409, 1952.

[52]
T.~M. Hospedales, A.~Antoniou, P.~Micaelli, and A.~J. Storkey.
\newblock Meta-learning in neural networks: A survey.
\newblock {\em IEEE transactions on pattern analysis and machine intelligence},
PP, 2021.

[53]
E.~Kobler, A.~Effland, K.~Kunisch, and T.~Pock.
\newblock Total deep variation for linear inverse problems.
\newblock In {\em Proceedings of the IEEE/CVF Conference on Computer Vision and
Pattern Recognition}, pages 7549--7558, 2020.

[54]
K.~Kunisch and T.~Pock.
\newblock A bilevel optimization approach for parameter learning in variational
\newblock {\em SIAM Journal on Imaging Sciences}, 6(2):938--983, 04 2013.

[55]
\newblock {\em Introduction to Smooth Manifolds}.
\newblock Graduate Texts in Mathematics. Springer, 2003.

[56]
A.~S. Lewis.
\newblock Active sets, nonsmoothness, and sensitivity.
\newblock {\em SIAM Journal on Optimization}, 13(3):702--725, 2002.

[57]
\newblock {\em Convergence rates of first-order operator splitting methods}.
\newblock PhD thesis, Normandie Universit{\'e}; GREYC CNRS UMR 6072, 2016.

[58]
J.~Liang, J.~Fadili, and G.~Peyré.
\newblock Activity identification and local linear convergence of
forward--backward-type methods.
\newblock {\em SIAM Journal on Optimization}, 27(1):408--437, 2017.

[59]
J.~Liang, J.~M. Fadili, and G.~Peyr\'{e}.
\newblock Local linear convergence of forward–backward under partial
\newblock In {\em Proceedings of the 27th International Conference on Neural
Information Processing Systems - Volume 2}, NIPS'14, pages 1970--1978,
Cambridge, MA, USA, 2014. MIT Press.

[60]
P.-L. Lions and B.~Mercier.
\newblock Splitting algorithms for the sum of two nonlinear operators.
\newblock {\em SIAM Journal on Numerical Analysis}, 16(6):964--979, 1979.

[61]
T.-T. Lu and S.-H. Shiou.
\newblock Inverses of 2$\times$ 2 block matrices.
\newblock {\em Computers \& Mathematics with Applications}, 43(1-2):119--129,

[62]
D.~Maclaurin, D.~Duvenaud, and R.~Adams.
\newblock Gradient-based hyperparameter optimization through reversible
\newblock In {\em International conference on machine learning}, pages
2113--2122. PMLR, 2015.

[63]
D.~Martin, C.~Fowlkes, D.~Tal, and J.~Malik.
\newblock A database of human segmented natural images and its application to
evaluating segmentation algorithms and measuring ecological statistics.
\newblock In {\em Proceedings of Eighth IEEE International Conference on
Computer Vision}, volume~2, pages 416--423. IEEE, 2001.

[64]
S.~Mehmood and P.~Ochs.
\newblock Automatic differentiation of some first-order methods in parametric
\newblock In S.~Chiappa and R.~Calandra, editors, {\em The 23rd International
Conference on Artificial Intelligence and Statistics, {AISTATS} 2020, 26-28
August 2020, Online [Palermo, Sicily, Italy]}, volume 108 of {\em Proceedings
of Machine Learning Research}, pages 1584--1594. PMLR, 2020.

[65]
A.~Mensch and M.~Blondel.
\newblock Differentiable dynamic programming for structured prediction and
\newblock In {\em International Conference on Machine Learning}, pages
3462--3471. PMLR, 2018.

[66]
\newblock On bayesian methods for seeking the extremum and their application.
\newblock In {\em IFIP Congress}, pages 195--200, 1977.

[67]
\newblock A method for unconstrained convex minimization problem with the rate
of convergence o (1/k\^{} 2).
\newblock In {\em Doklady an ussr}, volume 269, pages 543--547, 1983.

[68]
P.~Ochs, T.~Brox, and T.~Pock.
\newblock ipiasco: Inertial proximal algorithm for strongly convex
\newblock {\em Journal of Mathematical Imaging and Vision}, 53(2):171--181,
October 2015.

[69]
P.~Ochs, R.~Ranftl, T.~Brox, and T.~Pock.
\newblock Bilevel optimization with nonsmooth lower level problems.
\newblock In J.-F. Aujol, M.~Nikolova, and N.~Papadakis, editors, {\em
International Conference on Scale Space and Variational Methods in Computer
Vision (SSVM)}, volume 9087 of {\em Lecture Notes in Computer Science}, pages
654--665. Springer, 2015.

[70]
P.~Ochs, R.~Ranftl, T.~Brox, and T.~Pock.
\newblock Techniques for gradient based bilevel optimization with nonsmooth
lower level problems.
\newblock {\em Journal of Mathematical Imaging and Vision}, 56(2):175--194,

[71]
A.~Paszke, S.~Gross, F.~Massa, A.~Lerer, J.~Bradbury, G.~Chanan, T.~Killeen,
Z.~Lin, N.~Gimelshein, L.~Antiga, A.~Desmaison, A.~Kopf, E.~Yang, Z.~DeVito,
M.~Raison, A.~Tejani, S.~Chilamkurthy, B.~Steiner, L.~Fang, J.~Bai, and
\newblock Pytorch: An imperative style, high-performance deep learning library.
\newblock In H.~Wallach, H.~Larochelle, A.~Beygelzimer, F.~d\textquotesingle
Alch\'{e}-Buc, E.~Fox, and R.~Garnett, editors, {\em Advances in Neural
Information Processing Systems 32}, pages 8024--8035. Curran Associates,
Inc., 2019.

[72]
\newblock Hyperparameter optimization with approximate gradient.
\newblock In {\em International conference on machine learning}, pages
737--746. PMLR, 2016.

[73]
B.~Recht, M.~Fazel, and P.~A. Parrilo.
\newblock Guaranteed minimum-rank solutions of linear matrix equations via
nuclear norm minimization.
\newblock {\em SIAM review}, 52(3):471--501, 2010.

[74]
E.~S. Riis.
\newblock {\em Geometric numerical integration for optimisation.}
\newblock PhD thesis, University of Cambridge, 2020.

[75]
R.~Rockafellar and R.~J.-B. Wets.
\newblock {\em Variational Analysis}.
\newblock Springer Verlag, Heidelberg, Berlin, New York, 1998.

[76]
S.~Roth and M.~J. Black.
\newblock Fields of experts: A framework for learning image priors.
\newblock In {\em 2005 IEEE Computer Society Conference on Computer Vision and
Pattern Recognition (CVPR'05)}, volume~2, pages 860--867. IEEE, 2005.

[77]
L.~I. Rudin, S.~Osher, and E.~Fatemi.
\newblock Nonlinear total variation based noise removal algorithms.
\newblock {\em Physica D: Nonlinear Phenomena}, 60(1-4):259--268, 1992.

[78]
D.~E. Rumelhart, G.~E. Hinton, and R.~J. Williams.
\newblock Learning representations by back-propagating errors.
\newblock {\em Nature}, 323(6088):533--536, 1986.

[79]
S.~Schlenkrich, A.~Walther, N.~R. Gauger, and R.~Heinrich.
\newblock Differentiating fixed point iterations with {ADOL-C}: Gradient
calculation for fluid dynamics.
\newblock In {\em Modeling, Simulation and Optimization of Complex Processes},
pages 499--508. Springer, 2008.

[80]
\newblock Regression shrinkage and selection via the lasso.
\newblock {\em Journal of the Royal Statistical Society: Series B
(Methodological)}, 58(1):267--288, 1996.

[81]
S.~Vaiter, C.~Deledalle, J.~Fadili, G.~Peyr{\'e}, and C.~Dossal.
\newblock The degrees of freedom of partly smooth regularizers.
\newblock {\em Annals of the Institute of Statistical Mathematics},
69(4):791--832, 2017.

[82]
Y.~M. Volin and G.~Ostrovskii.
\newblock Automatic computation of derivatives with the use of the multilevel
differentiating technique—1. algorithmic basis.
\newblock {\em Computers \& mathematics with applications}, 11(11):1099--1114,

[83]
R.~J. Williams and J.~Peng.
\newblock An efficient gradient-based algorithm for on-line training of
recurrent network trajectories.
\newblock {\em Neural computation}, 2(4):490--501, 1990.

[84]
E.~Winston and J.~Z. Kolter.
\newblock Monotone operator equilibrium networks.
\newblock {\em Advances in neural information processing systems},
33:10718--10728, 2020.

[85]
M.~Yuan and Y.~Lin.
\newblock Model selection and estimation in regression with grouped variables.
\newblock {\em Journal of the Royal Statistical Society: Series B (Statistical
Methodology)}, 68(1):49--67, 2006.

\end{thebibliography}

\end{document}

%%% Local Variables:
%%% mode: latex
%%% TeX-master: t
%%% End:


# Experimental Test of Contextuality based on State Discrimination with a
Single Qubit

Qiuxin Zhang  Department of Physics, Renmin University of China, Beijing
100872, China  Chenhao Zhu  Department of Physics, Renmin University of China,
Beijing 100872, China  Yuxin Wang  Department of Physics, Renmin University of
China, Beijing 100872, China  Liangyu Ding  Department of Physics, Renmin
University of China, Beijing 100872, China  Tingting Shi  Department of
Physics, Renmin University of China, Beijing 100872, China  Xiang Zhang
Department of Physics, Renmin University of China, Beijing 100872, China
Beijing Academy of Quantum Information Sciences, Beijing 100193, China
Shuaining Zhang<EMAIL_ADDRESS>Department of Physics, Renmin
University of China, Beijing 100872, China  Beijing Academy of Quantum
Information Sciences, Beijing 100193, China  Wei Zhang<EMAIL_ADDRESS>Department of Physics, Renmin University of China, Beijing 100872, China
Beijing Academy of Quantum Information Sciences, Beijing 100193, China
Beijing Key Laboratory of Opto-electronic Functional Materials and Micro-nano
Devices, Renmin University of China, Beijing 100872, China

###### Abstract

Exploring quantum phenomena beyond predictions of any classical model has
fundamental importance to understand the boundary of classical and quantum
descriptions of nature. As a typical property that a quantum system behaves
distinctively from a classical counterpart, contextuality has been studied
extensively and verified experimentally in systems composed of at least three
levels (qutrit). Here we extend the scope of experimental test of
contextuality to a minimal quantum system of only two states (qubit) by
implementing the minimum error state discrimination on a single 171Yb+ ion. We
observe a substantial violation of a no-go inequality derived by assuming non-
contextuality, and firmly conclude that the measured results of state
discrimination cannot be reconciled with any non-contextual description. We
also quantify the contextual advantage of state discrimination and the
tolerance against quantum noises.

###### pacs:

03.65.Ta, 42.50.Dv, 03.65.Aa

Introduction.– In 1935, Einstein, Podolsky, and Rosen EPR35  firstly posed the
question of whether quantum observables have physical reality. It is then
realized that such a question has significant physical indication rather than
a conceptual or metaphysical notion. A variety of statements called no-go
theorems have been formulated in mathematical terms to exclude from our
natural world some different assumptions of reality. The pioneering works by
Bell Bell65 ; BrunnerRMP  and Kochen and Specker (KS) Kochen ; contextreview
provide two successful quantitative descriptions, where the violation of Bell
inequality effectively precludes local causation, and the no-go theorem of KS
theory proves that quantum mechanics will give conflicting predictions from
any classical non-contextual model Kochen ; contextreview . These statements
are of great importance since they are both rigorous in mathematics and
realistic in experiments PRA4 , where the test of Bell inequality usually
requires entangled states, and KS inequality can be violated for a more
general state of a quantum system.

The test of KS inequality has been demonstrated in various physical platforms
composed of photons photons1 ; photons2 , superconducting systems
NCSuperconduct , nuclear spins NVCenter , ions ions1 ; ions3 ; context1 , or
neutrons neutrons1 . In fact, in the conventional formalism of KS theory, a
qutrit system including at least three levels is the simplest case to show
quantum contextuality. Previous works are then demonstrated either for two-
qubit systems or a qutrit with multiple observables NVCenter ; ions1 ; Mermin
; Peres ; Cabello ; contextnature ; ions3 ; NCSuperconduct ; photons1 ;
neutrons1 ; photons2 ; kcbs ; Cabellokcbs ; context1 . Recently, a research
about contextual advantage for quantum state discrimination (SD) suggests that
a test of contextuality may be extended to even simpler systems and with fewer
observables PRX .

In this work, we show the first experimental test of quantum contextuality in
the simplest possible system of a single qubit with only two observable
measurements based on the non-classicality of state discrimination. State
discrimination focuses on how to effectively determine the state of a quantum
system according to a set of measured probabilities with a known prior
distribution. For a Hermitian system, one can prove that two non-orthogonal
pure states cannot be entirely discriminated within either a classical or a
quantum theory. However, a quantum description can significantly improves the
probability of success, such that the research along this line is of great
interest in both the study of fundamental quantum physics PBR  and the search
of potential applications in quantum information processing QKD ; QIP . This
advantage occurs exclusively in non-classical theory and has been studied in
various areas such as quantum clone theorem NJP1 ; NJP2 ; NJP6 , dimension
witness NJP18 , and mutual information NJP19 . Using a trapped 171Yb+ ion, we
realize for the first time an experimental test of contextuality under the
minimal error state discrimination (MESD) scenario with equal prior
probability, where a no-go inequality that tolerates experimental error in
non-contextual notation can be derived PRX . With high-fidelity single-qubit
logic gates and highly distinguishable state-dependent fluorescence detection,
a significant violation of the inequality is observed and the experimental
results are well agreed with theoretical predictions of a quantum description.

Preliminary.– We consider the problem of distinguishing two pure states
$\left|\psi\right\rangle$ and $\left|\phi\right\rangle$ with equal priori
probability, as shown in Fig. 1(a). The states $\left|\bar{\psi}\right\rangle$
and $\left|\bar{\phi}\right\rangle$ are orthogonal to
$\left|\psi\right\rangle$ and $\left|\phi\right\rangle$, respectively, and are
introduced as complement to construct operational equivalence PRA4  (See Data
analysis and results) and derive the no-go inequality with the hypothesis of
non-contextuality. By defining the Helstrom measurement Helstrom  basis set
(the dashed line) as
$\\{M_{d\left|\psi\right\rangle}=1/2(I+\sigma_{z}),M_{d\left|\phi\right\rangle}=1/2(I-\sigma_{z})\\}$
with $M_{d{\left|\psi\right\rangle}}+M_{d\left|\phi\right\rangle}=I$, we can
write the probability of a successful differentiation between
$\left|\psi\right\rangle$ and $\left|\phi\right\rangle$ as $s=\frac{1}{2}({\rm
Tr}[M_{d\left|\psi\right\rangle}\left|\psi\right\rangle\left\langle\psi\right|]+{\rm
Tr}[M_{d\left|\phi\right\rangle}\left|\phi\right\rangle\left\langle\phi\right|])$.
It is reasonable and inevitable for maximizing the successful probability to
assume that the measurements have the symmetry property ${\rm
Tr}[M_{d\left|\psi\right\rangle}\left|\psi\right\rangle\left\langle\psi\right|]={\rm
Tr}[M_{d\left|\phi\right\rangle}\left|\phi\right\rangle\left\langle\phi\right|]$,
so that $s={\rm
Tr}[M_{d\left|\psi\right\rangle}\left|\psi\right\rangle\left\langle\psi\right|]={\rm
Tr}[M_{d\left|\phi\right\rangle}\left|\phi\right\rangle\left\langle\phi\right|]$.
With that, the quantity $s$ is related to the overlap of
$\left|\phi\right\rangle$ and $\left|\psi\right\rangle$ as

$\displaystyle
s=\frac{1}{2}\left(1+\sqrt{1-|\langle\phi|\psi\rangle|^{2}}\right)=\frac{1}{2}\left(1+\sqrt{1-c}\right),$
(1)

where the overlap $c\equiv|\langle\phi|\psi\rangle|^{2}$ is usually referred
as confusability. It is obvious that when the direction of the measurement
basis $M_{d\left|\psi(\phi)\right\rangle}$ and the direction of the vector
$(\left|\psi\right\rangle-\left|\phi\right\rangle)$ are parallel as shown in
Fig. 1(a), the error rate of the two states discrimination is the lowest, by
which we call the minimum-error state discrimination (MESD).

Figure 1:  (a) The scheme of state preparations on the $x$-$z$ Bloch plane.
The state $\left|\psi(\phi)\right\rangle$ (dark arrow) has an angle $\theta$
with measurement basis $M_{d\left|\psi(\phi)\right\rangle}$ (dashed line). We
change the value of $\theta$ from 0 to $\pi/2$ to prepare different state
ensembles, such that the two states change from orthogonal to parallel. The
state $\left|\bar{\psi}(\bar{\phi)}\right\rangle$ is orthogonal to the state
$\left|\psi(\phi)\right\rangle$. (b) The measurement scheme for states on the
$x$-$z$ Bloch plane. The ideal measurement basis are shown with gray dashed
lines. A measurement error $\epsilon$ is considered in the basis
$M^{\alpha}_{\left|\psi(\phi)\right\rangle}$($\alpha\in[0,\pi/2]$), which has
an angle deviation $\alpha$ from the ideal basis
$M_{\left|\psi(\phi)\right\rangle}$. The measurement basis
$M^{\alpha}_{\left|\phi(\psi)\right\rangle}$ and
$M^{\alpha}_{{\left|\bar{\phi}(\bar{\psi})\right\rangle}}$ are orthogonal. (c)
The trapped single 171Yb+ ion system with a four-rod Paul trap. The purple
beam is the 369.5nm laser for cooling, state initialization and state readout.
The single-qubit gate is realized by a microwave horn antenna. (d) The energy
level diagram of 171Yb+.

As proved in Ref. PRX , any ontological model of an operational theory with
non-contextuality must lead to the following no-go theorem for MESD

$\displaystyle s\leq 1-\frac{c}{2}.$ (2)

This inequality would be relaxed owing to the error of quantum operations or
measurement basis under realistic experimental conditions. To this end, one
obtains

$\displaystyle s\leq 1-\frac{c-\epsilon}{2},$ (3)

where $\epsilon$ is a quantitative indicator of measurement incorrectness,
which is introduced just like the one appeared in the Clauser-Horne-Shimony-
Holt (CHSH) inequality CHSH  to describe the error of measurement in the Bell
inequality. In our experiment, $\epsilon$ is related to the deviation of the
actual measurement basis
$\\{M_{\left|\psi\right\rangle}^{\alpha},M_{\left|\phi\right\rangle}^{\alpha}\\}$
[solid lines in Fig. 1(b)] from the ideal measurement axes (dashed lines) with
an angle $\alpha\in[0,\pi/2]$. We further choose the directions of deviation
for $M_{\left|\psi\right\rangle}^{\alpha}$ and
$M_{\left|\phi\right\rangle}^{\alpha}$ are opposite to each other, i.e., they
are tilted either away or close to the $x$-axis simultaneously.

Taking the measurement error into account, the quantum bound of success
probability of SD can be derived as PRX

$\displaystyle
s=\frac{1}{2}\left[1+\sqrt{1-\epsilon+2\sqrt{\epsilon(1-\epsilon)c(1-c)}+c(2\epsilon-1)}\right].$
(4)

Apparently, the expression of Eq. (4) reduces to that of Eq. (1) when
$\epsilon=0$, and the no-go result of Eq. (2) is the corresponding limiting
case of Eq. (3). The success probability $s$ of discrimination in quantum
description is higher than the upper bound of non-contextual theories. We use
the same terminology as in Ref. PRX  to label the measurements and outcomes in
a natural way. For instance, the outcome of
$M_{\left|\phi\right\rangle}^{\alpha}$ is more likely to occur than the
outcome of $M_{\left|\psi\right\rangle}^{\alpha}$ if one begins with
$P_{\left|\phi\right\rangle}$ rather than $P_{\left|\bar{\psi}\right\rangle}$
preparation, and so on. Thus, Eqs. (3) and (4) must satisfy the constraint

$\displaystyle\epsilon\leq c\leq 1-\epsilon.$ (5)

We also calculate the theoretical predications of
$c=\langle\psi|M^{\alpha}_{\left|\phi\right\rangle}|\psi\rangle$,
$s=\langle\phi|M_{d\left|\phi\right\rangle}|\phi\rangle$, and
$1-\epsilon=\langle\phi|M^{\alpha}_{\left|\phi\right\rangle}|\phi\rangle$
under a quantum description for different state preparations and measurements.
The precise predictions are listed in Table 1. Note that two additional
preparations along the $y$-axis are introduced as an extra dimension to
correct the experimental preparation inaccuracy and supplement the measurement
observable $\sigma_{y}$ to ensure the completeness of tomographic measurement.

Table 1: Theoretical predictions for different state preparations and
measurements. Each pair of preparation $P_{i}^{\theta}$ (with
$i\in\\{\left|\phi\right\rangle,\left|\bar{\phi}\right\rangle,\left|\psi\right\rangle,\left|\bar{\psi}\\}\right\rangle$)
and measurements$M_{j}$ (with
$M_{j}\in\\{M_{\left|\psi\right\rangle}^{\alpha},M_{\left|\phi\right\rangle}^{\alpha},M_{d\left|\phi\right\rangle}\\}$)
yields a set of parameters $c$ = $\sin^{2}[({\alpha-2\theta})/{2}]$,
$\epsilon$ = $\sin^{2}({\alpha}/{2})$ and $s$ = $\cos^{2}({\theta}/{2})$. To
ensure that the preparations and measurements are tomographically complete, we
add the preparations {$P_{y}$,$P_{\bar{y}}$} and measurement $M_{y}$ in the
$y$-axis.

| $P_{\left|\psi\right\rangle}^{\theta}$ | $P_{\left|\bar{\psi}\right\rangle}^{\theta}$ | $P_{\left|\phi\right\rangle}^{\theta}$ | $P_{\left|\bar{\phi}\right\rangle}^{\theta}$ | $P_{y}$ | $P_{\bar{y}}$
---|---|---|---|---|---|---
$M_{\left|\psi\right\rangle}^{\alpha}$ | $1-\epsilon$ | $\epsilon$ | $c$ | $1-c$ | 0.5 | 0.5
$M_{\left|\phi\right\rangle}^{\alpha}$ | $c$ | $1-c$ | $1-\epsilon$ | $\epsilon$ | 0.5 | 0.5
$M_{d\left|\phi\right\rangle}$ | $1-s$ | $s$ | $s$ | $1-s$ | 0.5 | 0.5
$M_{y}$ | 0.5 | 0.5 | 0.5 | 0.5 | 1 | 0

To test the no-go inequality Eq. (3), we need to prepare all states and
measure all observables as given in Table 1, under the condition of
operational equivalence. Then we can obtain the values of $\\{c,\epsilon,s\\}$
from the relations given above and substitute into Eq. (3). If a violation of
such an inequality is observed beyond experimental error, one can conclude
that no non-contextual model can account for the experimental observation.
Obviously this situation will be more difficult to occur with increasing
noise. In this sense, the no-go inequality Eq. (3) provides a criterion for
the boundary between quantum and classical regimes.

Experimental setup.– The experiment is performed on a single ${}^{171}{\rm
Yb}^{+}$ ion trapped in a four-rod Paul trap, as shown in Fig. 1(c). A pair of
clock states $\left|F=0,m_{F}=0\right\rangle$ and
$\left|F=1,m_{F}=0\right\rangle$ of the 2S1/2 ground-state manifold are used
to encode a qubit as {$\left|\downarrow\right\rangle$,
$\left|\uparrow\right\rangle$} with transition frequency
$\omega_{0}=2\pi\times 12.64281$ GHz [Fig. 1(d)]. The qubit can be initialized
by Doppler cooling and optical pumping, and detected by shining a 369.5 nm
laser to excite the state $\left|\uparrow\right\rangle$ to the 2P1/2 manifold
and collecting the fluorescent photons with a photo-multiplier tube (PMT)
Olmschenk . The error of detecting $\left|\uparrow\right\rangle$ for
$\left|\downarrow\right\rangle$ is 1.71% and the other one is 2.08%. We
calibrate the results to remove the detection errors as proposed in Ref.
Detection .

The states $\left|\psi\right\rangle$ and $\left|\phi\right\rangle$ that we
want to prepare lie in the $x$-$z$ plane, and can be denoted as

$\displaystyle\begin{aligned}
\left|\psi\right\rangle&=\operatorname{\cos}\left({\theta}/{2}\right)\left|\downarrow\right\rangle+\operatorname{\sin}\left({\theta}/{2}\right)\left|\uparrow\right\rangle,\\\
\left|\phi\right\rangle&=\operatorname{\cos}\left[{(\pi-\theta)}/{2}\right]\left|\downarrow\right\rangle-\operatorname{\sin}\left[{(\pi-\theta)}/{2}\right]\left|\uparrow\right\rangle\end{aligned}$
(6)

on the basis of
$\\{\left|\downarrow\right\rangle,\left|\uparrow\right\rangle\\}$. When
$\theta$ changes from 0 to $\pi/2$, the two states become closer to the
$x$-axis on the Bloch sphere with increasing confusability $c$. To prepare
such states, single-qubit quantum gates are implemented via microwave with
resonant frequency. By controlling the duration and phase of the microwave, we
can achieve a native single-qubit rotation

$\displaystyle
R(\beta,\gamma)=\left(\begin{array}[]{cc}\cos(\beta/2)&-ie^{-i\gamma}\sin(\beta/2)\\\
-ie^{i\gamma}\sin(\beta/2)&\cos(\beta/2)\end{array}\right),$ (9)

where $\beta$ and $\gamma$ are the rotation angle and phase, respectively. The
$\pi$ time of the Rabi oscillation is 7.51 ${\mu}s$, corresponding to a Rabi
frequency of $2\pi\times 66.6$ kHz.

The experimental process is as follows: After $500\ \mu s$ Doppler cooling,
the qubit is initialized to the $\left|\downarrow\right\rangle$ state by a $6\
\mu s$ optical pumping. Then, the set of states
$\\{P_{\left|\psi\right\rangle}^{\theta},P_{\left|\bar{\psi}\right\rangle}^{\theta},P_{\left|\phi\right\rangle}^{\theta},P_{\left|\bar{\phi}\right\rangle}^{\theta}\\}$
are prepared through single-qubit rotations
$\\{R_{y}(\theta),R_{y}(\pi+\theta),R_{y}(\pi-\theta),R_{y}^{\dagger}(\theta)\\}$,
where $R_{y}(\theta)=R(\theta,\pi/2)$ and
$R_{y}^{\dagger}(\theta)=R(\theta,3\pi/2)$. The observable set is chosen as
$\\{M_{\left|\psi\right\rangle}^{\alpha},M_{\left|\phi\right\rangle}^{\alpha},M_{d\left|\phi\right\rangle}\\}$,
which are realized by adding single-qubit operations
$\\{R_{y}^{\dagger}(\theta-\alpha),R_{y}^{\dagger}(\pi-\theta+\alpha),R_{y}^{\dagger}(\pi)\\}$
to rotate the measurement basis from the original $z$-basis to the desired
one. By scanning the parameter region $\\{\theta,\alpha\\}$, we can obtain the
dataset $\\{c,\epsilon,s\\}$.

The parameter $\epsilon$ is a function of the variable $\alpha$ as $\epsilon$
= $\sin^{2}({\alpha}/{2})$. Similarly, we have $c$ =
$\sin^{2}[({\alpha-2\theta})/{2}]$. When $\epsilon$ has a certain value, $c$
must meet the inequality Eq. (5), which can be reduced to

$\displaystyle\sin^{2}\left({\alpha}/{2}\right)\leq\sin^{2}\left[({\alpha-2\theta})/{2}\right]\leq\cos^{2}\left({\alpha}/{2}\right).$
(10)

Note that $\sin^{2}[({\alpha-2\theta})/{2}]\leq\cos^{2}({\alpha}/{2})$ is
always true with $\theta,\alpha\in[0,\pi/2]$. Therefore, when $\alpha$ is
fixed, the range of $\theta$ can be determined by solving
$\sin^{2}({\alpha}/{2})\leq\sin^{2}[({\alpha-2\theta})/{2}]$ with
$\theta,\alpha\in[0,\pi/2]$. In our experiment, the parameters $\theta$ and
$\alpha$ are scanned with an interval of 0.1, resulting in 136 sets of data
$\\{c,\epsilon,s\\}$. We repeat the experiment for 1000 trials for each choice
of the data set.

To correct the state preparation inaccuracy on the Bloch sphere and ensure the
measurements are tomographically complete, we also supplement extra state
preparations and observable measurements through the similar quantum
operations. Details of this procedure will be given in the subsequent
discussion and Supplemental Materials supp . As a result, the total state
preparations are
$\\{P_{\left|\psi\right\rangle}^{\theta},P_{\left|\bar{\psi}\right\rangle}^{\theta},P_{\left|\phi\right\rangle}^{\theta},P_{\left|\bar{\phi}\right\rangle}^{\theta},P_{y},P_{\bar{y}}\\}$,
and observable measurements are
$\\{M_{\left|\psi\right\rangle}^{\alpha},M_{\left|\phi\right\rangle}^{\alpha},M_{d\left|{\phi}\right\rangle},M_{y}\\}$,
as shown in Table 1. We emphasize that although more than two state
preparations and observable measurements are carried in the experiment, only
two observables with binary outcome are obtained from the set
$\\{M_{\left|\psi\right\rangle}^{\alpha},M_{\left|\phi\right\rangle}^{\alpha},M_{d\left|{\phi}\right\rangle}\\}$,
and are subsequently used to test the no-go inequality and show the boundary
between quantum and classical regimes.

Figure 2: The results of data analysis for operational equivalence. (a) The
states
$\\{\left|\phi\right\rangle,\left|\bar{\phi}\right\rangle,\left|\psi\right\rangle,\left|\bar{\psi}\right\rangle\\}$
prepared in experiments are plotted on the $x$-$z$ plane of the Bloch sphere
with $\theta$ = 1.1 rad and $\alpha$ = 0 rad, corresponding to the point with
index 12 in panels (b), (c) and (d). Raw preparations (blue triangles) are
given from the experimental data. Primary preparations (orange dots) are the
results of fitting with the GPT theory. Secondary preparations (green
rectangles) are chosen to satisfy the assumption of operational equivalence.
(b), (c), and (d) show that the probabilities
$\mathcal{P}_{\frac{1}{2}\left|\psi\right\rangle+\frac{1}{2}\left|\bar{\psi}\right\rangle,M_{i}}$
reach statistical operational equivalence, i.e.
$\mathcal{P}_{\frac{1}{2}\left|\psi\right\rangle+\frac{1}{2}\left|\bar{\psi}\right\rangle,M_{i}}=\frac{1}{2}$
after data processing. The results of
$\mathcal{P}_{\frac{1}{2}\left|\phi\right\rangle+\frac{1}{2}\bar{\left|\phi\right\rangle},M_{i}}$
are given in the Supplemental Materials.

Data analysis and results.– The test of no-go inequality derived from the non-
contextual model requires the data meet operational equivalence in a general
framework regardless of quantum or classic. For this purpose, we need to
process the raw experimental data in two stages NC ; Mikeprxq  (See details in
the Supplemental Materials). The first phase fits the data of raw preparations
using generalized probability theory (GPT) to ensure the states are not
exclusively dependent on quantum theory under the condition of complete
tomographic measurements. The results of this stage, refereed as the primary
preparations, are then supplied to the second phase which treats the data to
perfectly comply to operational equivalence, i.e.

$\displaystyle\mathcal{P}_{\frac{1}{2}\phi+\frac{1}{2}\bar{\phi},M_{i}}=\mathcal{P}_{\frac{1}{2}\left|\psi\right\rangle+\frac{1}{2}\left|\bar{\psi}\right\rangle,M_{i}}=\frac{1}{2}.$
(11)

Here,
$\mathcal{P}_{\frac{1}{2}\left|\phi\right\rangle+\frac{1}{2}\left|\bar{\phi}\right\rangle,M_{i}}$
denotes a process associated with a particular
$M_{i}\in\\{M_{\left|\phi\right\rangle},M_{\left|\psi\right\rangle},M_{d\left|\phi\right\rangle}\\}$,
which is randomly selected from the average random mid-sample of the
$P_{\left|\phi\right\rangle}$ and $P_{\left|\bar{\phi}\right\rangle}$
preparation process in theory. The case is analogous for
$\mathcal{P}_{\frac{1}{2}\psi+\frac{1}{2}\bar{\psi},M_{i}}$. In experiment,
$\mathcal{P}_{\frac{1}{2}\left|\phi\right\rangle+\frac{1}{2}\left|\bar{\phi}\right\rangle,M_{i}}$
is realized by a linear combination of the statistics of
$\left|\phi\right\rangle$ and $\left|\bar{\phi}\right\rangle$ under
measurement $M_{i}$, i.e.,

$\displaystyle\mathcal{P}_{\frac{1}{2}\left|\phi\right\rangle+\frac{1}{2}\left|\bar{\phi}\right\rangle,M_{i}}=\frac{1}{2}\left[P\left(1\mid
P_{\left|\phi\right\rangle},M_{i}\right)+P\left(1\mid
P_{\left|\bar{\phi}\right\rangle},M_{i}\right)\right].$ (12)

Here, $P(1|P_{\phi(\bar{\phi})},M_{i})$ denotes the frequency when the
measurement axis is $M_{i}$ and the outcome is 1. Using maximum likelihood
estimation method, we obtain the secondary state preparations by minimizing
the target function as the summation of distances from their corresponding
primaries.

The two stages of data processing alter the position of states in the $x$-$z$
plane as illustrated in Fig. 2(a). The data processing outcomes are plotted in
Figs. 2(b)-2(d). Notice that the raw and primary data of state preparations
are slightly off 1/2 owing to experimental errors in operation and
measurement. The primary preparations fitted from the GPT theory have a
slightly larger variance, which can be hardly observed in Figs. 2(b)-2(d). The
results of secondary preparations show that all the 136 points fulfill the
operational equivalence criterion with an exactly equal distribution 1/2. The
similarity between states before and after data processing is evaluated as the
average of mixture weights, which is determined as $97.13\%\pm 2.3\%$ in our
experiment. Further details and results can be found in the Supplemental
Materials.

Figure 3: The results of contextual advantage for state discrimination. (a)
The results of success probability $s$ versus confusability $c$ for
$\epsilon$= $\\{0,0.022,0.146,0.315\\}$. The points of dots, squares, diamonds
and triangles are the experiment data. The dashed lines represent the non-
contextual upper bound $s_{\rm non-con}$ given by Eq. (3), and the solid lines
are the predictions by a quantum theory $s_{\rm theory}$. The shaded areas
mark the contextual advantage $\Delta s_{\rm theory}=s_{\rm theory}-s_{\rm
non-con}$ for state discrimination. (b) The results of $\Delta s$ versus $c$
with the same $\epsilon$’s and line shapes as in (a). The red pentagram and
hexagonal star are the highest points with and without $\epsilon$, locating at
($c,\epsilon$)=(0.853, 0.146) and ($c$,$\epsilon$)=(0.75,0), respectively. (c)
and (d) The false-color plots of $\Delta s$ by changing $c$ and $\epsilon$.

We use the secondary preparation data which satisfies operational equivalence
to calculate the results of $\\{c,\epsilon,s\\}$ and show the violation of the
no-go inequality Eq. (3). In Fig. 3(a), we plot the probability of successful
state discrimination $s$ by varying the confusability $c$ with fixed
measurement error
$\epsilon=\\{\sin(0/2)^{2},\sin(0.3/2)^{2},\sin(0.785/2)^{2},\sin(1.2/2)^{2}\\}=\\{0,0.022,0.146,0.315\\}$.
The measured $s_{\rm exp}$ is clearly above the upper bound $s_{\rm non-
con}=1-(c-\epsilon)/2$ for any non-contextual model (dashed lines), and agree
well with theoretical predictions $s_{\rm theory}$ of Eq. (4) given by quantum
theory with contextuality (solid lines). The shaded region between the dashed
and solid lines thus denotes the advantage of a contextual theory in SD. We
also notice that in the limiting case of $c=\epsilon$, the success probability
of SD can reach $s=1$ even for a pair non-orthogonal state preparations with
$c>0$. This is because in the presence of a finite $\epsilon$, a small
confusability between two preparations with $c\leq\epsilon$ can always be
attributed to the measurement error, such that the two states are perfectly
distinguished within an ontological description of the system, regardless of
contextuality or non-contextuality.

To show the contextual advantage of SD more clearly, we plot in Fig. 3(b) the
exceeding of $s$ over the non-contextual bound $\Delta s_{\rm
theory/exp}=s_{\rm theory/exp}-s_{\rm non-con}$ for each choice of $\epsilon$.
We emphasize that the highest value of $\Delta s_{\rm theory}=s_{\rm
theory}-s_{\rm non-con}$ is 0.207 for the parameters ($c,\epsilon$)=(0.853,
0.146), which clearly exceeds the maximal value of 0.125 for the case with no
measurement error achieved at ($c$,$\epsilon$)=(0.75,0). These two points are
marked in Fig. 3(b) with a red pentagram and a hexagonal star, respectively.
It is obvious that the addition of a finite $\epsilon$ significantly improves
the contextual advantage of SD over any non-contextual model. In experiment,
the location of the pentagram and hexagonal star are determined at
$\\{c,\epsilon,s\\}=\\{0.750\pm 0.010,0.0000\pm 0.0006,0.125\pm 0.007\\}$ and
$\\{c,\epsilon,s\\}=\\{0.853\pm 0.009,0.146\pm 0.011,0.206\pm 0.006\\}$,
respectively. The quantity $\Delta s$ also manifests the tolerance against
noise in SD. For instance, one can prove that in a contextual theory $s$
decreases with increasing depolarization error of state preparation and
measurement PRX , which eventually compromises the violation of inequality Eq.
3. Thus, a higher value of $\Delta s$ achieved at a finite $\epsilon$
indicates that a quantum SD is more robust against depolarization noises.

The results of $\Delta s$ for all the 136 points obtained by theory and
experiment are plotted in false color in Figs. 3(c) and 3(d), respectively. To
quantify the resemblance between theory and experiment, we define the average
fidelity as

$\displaystyle f=1-\frac{1}{N}\sum_{i=1}^{N=136}\frac{\left|\Delta s_{{i,\rm
exp}}-\Delta s_{i,{\rm theory}}\right|}{\Delta s_{i,{\rm theory}}},$ (13)

and obtains $f=98.82\%\pm 0.12\%$. Since the gate infidelity induced by the
fluctuation of microwave parameters is very low, which is $\sim(2.0\pm
0.3)\times$10-4 as characterized by randomized benchmarking method RB , the
main error source is the sampling error with an average of $1.29\%\pm 0.45\%$.
(See details in the Supplemental Materials).

Conclusion.– We demonstrate an experimental no-go test of contextuality
derived within the frame work of minimum error state discrimination in a
single qubit system with only two observable measurements. With high-
fidelities of state initialization, manipulation and detection, we observe a
probability of successful differentiation between two non-orthogonal state
preparations which significantly exceeds the upper bound for any non-
contextual theories. Our results represent the first experimental ruling out
of non-contextual descriptions in a simple qubit system with only two
observables. We also introduce the measurement error to quantify the
contextual advantage and the tolerance over noise of quantum state
discrimination. Our work not only deepen the understanding of the quantum
nature of contextuality, but also motivate future research on noise effect
from quantitative aspects.

###### Acknowledgements.

We are thankful for support from the Beijing Natural Science Foundation (Grant
No. Z180013), the National Natural Science Foundation of China (Grants No.
12074427, and No. 12074428), and the National Key R&D Program of China (Grant
No. 2018YFA0306501). S. Zhang also thanks the China Postdoctoral Science
Foundation (Grant No. BX20200379 and No. 2021M693478) for support.

## References

  * (1) A. Einstein, B. Podolsky, and N. Rosen, Can quantum-mechanical description of physical reality be considered complete? Phys. Rev. 47, 777 (1935).
  * (2) J. S. Bell, On the Einstein Podolsky Rosen paradox, Physics 1, 195 (1964).
  * (3) N. Brunner, D. Cavalcanti, S. Pironio, V. Scarani, and S. Wehner, Bell nonlocality, Rev. Mod. Phys. 86, 419 (2014).
  * (4) S. Kochen and E. P. Specker, The Problem of Hidden Variables in Quantum Mechanics, J. Math. Mech. 17, 59 (1967).
  * (5) C. Budroni, A. Cabello, O. Güuhne, M. Kleinmann, and J.-Å. Larsson, Quantum contextuality, arXiv:2102.13036.
  * (6) W. Spekkens, Contextuality for Preparations, Transformations, and Unsharp Measurements, Phys. Rev. A 71, 052108 (2005).
  * (7) A. Zhang, H. Xu, J. Xie, H. Zhang, B. J. Smith, M. S. Kim, and L. Zhang, Experimental test of contextuality in quantum and classical systems, Phys. Rev. Lett. 122, 080401 (2019).
  * (8) D. Qu, K. Wang, L. Xiao, X. Zhan, and P. Xue, State-independent test of quantum contextuality with either single photons or coherent light, NPJ Quantum Inf. 7, 1 (2021).
  * (9) M. Jerger, Y. Reshitnyk, M. Oppliger, A. Potočnik, M. Mondal, A. Wallraff, K. Goodenough, S. Wehner, K. Juliusson, N. K. Langford, and A. Fedorov, Contextuality without nonlocality in a superconducting quantum system, Nat. Commun. 7, 1 (2016).
  * (10) S. B. van Dam, J. Cramer, T. H. Taminiau, and R. Hanso, Multipartite entanglement generation and contextuality tests using nondestructive three-qubit parity measurements, Phys. Rev. Lett. 123, 050401 (2019).
  * (11) F. M. Leupold, M. Malinowski, C. Zhang, V. Negnevitsky, A. Cabello, J. Alonso, and J. P. Home, Sustained state-independent quantum contextual correlations from a single ion, Phys. Rev. Lett. 120, 180401 (2018).
  * (12) P. F. Wang, J. H. Zhang, C. Y. Luan, Mark Um, Y. Wang, M. Qiao, T. Xie, J. N. Zhang, A. Cabello, and K. Kim, Significant loophole-free test of Kochen-Specker contextuality using two species of atomic ions, Sci. Adv. 8, eabk1660 (2022).
  * (13) X. Zhang, M. Um, J. Zhang, S. An, Y. Wang, D.-l. Deng, C. Shen, L.-M. Duan, and K. Kim, State-independent experimental test of quantum contextuality with a single trapped ion, Phys. Rev. Lett. 110, 070401 (2013).
  * (14) Y. Hasegawa, R. Loidl, G. Badurek, M. Baron, and H. Rauch, Violation of a Bell-like inequality in single-neutron interferometry, Nature 425, 45 (2003).
  * (15) N. D. Mermin, Extreme quantum entanglement in a superposition of macroscopically distinct states, Phys. Rev. Lett. 65, 1838 (1990).
  * (16) A. Peres, Two simple proofs of the Kochen-Specker theorem, J. Phys. A 24, L175(1991).
  * (17) A. Cabello, Experimentally Testable State-Independent Quantum Contextuality, Phys. Rev. Lett. 101, 210401 (2008).
  * (18) G. Kirchmair, F. Zähringer, R. Gerritsma, M. Kleinmann, O. Gühne, A. Cabello, R. Blatt, and C. F. Roos, State-independent experimental test of quantum contextuality, Nature 460, 494 (2009).
  * (19) A. A. Klyachko, M. A. Can, S. Binicioǧlu, and A. S. Shumovsky, Simple test for hidden variables in spin-1 systems, Phys. Rev. Lett. 101, 020403 (2008).
  * (20) M. Kleinmann, C. Budroni, J-Å. Larsson, O. Gühne, and A. Cabello, Optimal inequalities for state-independent contextuality, Phys. Rev. Lett. 109, 250402 (2012).
  * (21) D. Schmid and R.W. Spekkens, Contextual Advantage for State Discrimination, Phys. Rev. X 8, 011015 (2018).
  * (22) M. F. Pusey, B. Jonathan, and R. Terry, On the reality of the quantum state, Nat. Phys. 8, 475 (2012).
  * (23) C. H. Bennett, Quantum cryptography using any two nonorthogonal states, Phys. Rev. Lett. 68, 3121 (1992).
  * (24) J. Bae and L. C. Kwek, Quantum state discrimination and its applications, J. Phys. A 48, 083001 (2015).
  * (25) C. H. Bennett and G. Brassard, Quantum cryptography, Proc. IEEE Int. Conf. on Computers, Systems and Signal Processing, Bangalore, India. 175 (1984).
  * (26) W. K. Wootters and W. H. Zurek, A single quantum cannot be cloned, Nature 299, 802 (1982).
  * (27) N. Gisin, Quantum cloning without signaling, Phys. Lett. A 242, 1 (1998).
  * (28) N. Brunner, M. Navascués, and T. Vértesi, Dimension witnesses and quantum state discrimination, Phys. Rev. Lett. 110, 150501 (2013).
  * (29) R. Koenig, R. Renner, and C. Schaffner, The operational meaning of min-and max-entropy, IEEE Trans. Inf. Theory 55, 4337 (2009).
  * (30) C. W. Helstrom, Quantum Detection and Estimation Theory, J. Stat. Phys. 1, 231 (1969).
  * (31) J. F. Clauser, M. A. Horne, A. Shimony, and R. A. Holt, Proposed Experiment to Test Local Hidden-Variable Theories, Phys. Rev. Lett. 23, 880 (1969).
  * (32) S. Olmschenk, K. C. Younge, D. L. Moehring, D. N. Matsukevich, P. Maunz, and C. Monroe, Manipulation and detection of a trapped Yb hyperfine qubit, Phys. Rev. A 76, 052314 (2007).
  * (33) C. Shen and L. M. Duan, Correcting detection errors in quantum state engineering through data processing, New J. Phys. 14, 053053 (2012).
  * (34) See Supplemental Materials for the data and error analysis.
  * (35) M. D. Mazurek, M. F. Pusey, R. Kunjwal, K. J. Resch, and R. W. Spekkens, An experimental test of noncontextuality without unphysical idealizations, Nat. Commun. 7, 11780 (2016).
  * (36) M. D. Mazurek, M. F. Pusey, K. J. Resch, and R. W. Spekkens, Experimentally bounding deviations from quantum theory in the landscape of generalized probabilistic theories, PRX Quantum 2, 020302 (2021).
  * (37) E. Knill, D. Leibfried, R. Reichle, J. Britton, R. B. Blakestad, J. D. Jost, C. Langer, R. Ozeri, S. Seidelin, and D. J. Wineland, Randomized benchmarking of quantum gates, Phys. Rev. A 77, 012307 (2008).