# Chance-Constrained Optimization in Contact-Rich Systems
for Robust Manipulation

Yuki Shirai†, Devesh K. Jha‡, Arvind U. Raghunathan‡, Diego Romeres‡ † Yuki
Shirai is with the Department of Mechanical and Aerospace Engineering,
University of California, Los Angeles, CA, USA 90095
<EMAIL_ADDRESS>K. Jha, Arvind U. Raghunathan and Diego
Romeres are with Mitsubishi Electric Research Laboratories (MERL), Cambridge,
MA, USA 02139<EMAIL_ADDRESS>
###### Abstract

This paper presents a chance-constrained formulation for robust trajectory
optimization during manipulation. In particular, we present a chance-
constrained optimization for Stochastic Discrete-time Linear Complementarity
Systems (SDLCS). To solve the optimization problem, we formulate Mixed-Integer
Quadratic Programming with Chance Constraints (MIQPCC). In our formulation, we
explicitly consider joint chance constraints for complementarity as well as
states to capture the stochastic evolution of dynamics. We evaluate robustness
of our optimized trajectories in simulation on several systems. The proposed
approach outperforms some recent approaches for robust trajectory optimization
for SDLCS.

##  I Introduction

Contacts are central to manipulation problems. Consequently, contact modeling
has been an active area of research in robotics since the last several decades
[1, 2, 3, 4, 5]. One of the most popular approaches to model contact dynamics
is using Linear Complementarity Problem (LCP). LCP models are widely used for
modeling contact dynamics in academia as well as in several physics simulation
engines such as Bullet, ODE, etc. Trajectory optimization (TO) of LCP-based
contact models has been used for manipulation [6, 7] and legged locomotion
[8]. Lyapunov stability of linear systems with complementarity systems has
also been studied [9, 10, 11]. However, most of these works assume
deterministic contact models to perform TO. In reality, frictional interaction
systems suffer from several uncertainties which lead to stochastic dynamics
and thus, it is important to consider uncertainty during TO. Modeling
uncertainty in LCP-based contact models leads to Stochastic Discrete-time
Linear Complementarity System (SDLCS).

Figure 1: This paper presents chance-constrained optimization for SDLCS. The
figure shows the case of a sliding box on a plane where the coefficient of
friction is a Gaussian random variable. Note that $w$ and $v$ are additive
uncertainty terms.

Fig. 1 pictorially shows a SDLCS we study in this paper. We consider the SDLCS
that has uncertainty in parameters and additive noises in dynamics and
complementarity constraints. As shown in Fig. 1, one should notice that
uncertainty leads to stochastic evolution of system states in SDLCS. Thus, a
robust optimization formulation should consider the uncertainty in state
evolution. In some recent works that consider stochastic complementarity
constraints, an expected residual minimization (ERM)-based [12] penalty is
used to solve the robust optimization problem [13]. A major shortcoming of
such an approach is that it fails to capture the stochastic evolution of
system dynamics due to the stochastic complementarity constraint. In [14], the
authors augment the formulation in [13] with chance constraints. However, this
formulation has certain fundamental shortcomings which prevent constraint
satisfaction guarantees. We present a formulation that circumvents these
shortcomings by using a mixed integer formulation. Using some relaxation of
the original joint chance complementarity constraint problem for the SDLCS,
the resulting problem can be solved using mixed integer programming.

In this paper, we present a formulation of robust trajectory optimization for
SDLCS. Since worst-case robust optimization is quite conservative and does not
explicitly discuss stochastic evolution of states [15], this work considers
probabilistic optimization with stochastic evolution of states. Robustness to
uncertainty is provided by enforcing probabilistic satisfaction of state
constraints. Under certain simplifications, we show that the chance-
constrained problem can be reformulated as a Mixed Integer Quadratic Program
with Chance Constraints.

Contributions. This paper has the following contributions:

  1. 1.

We present a novel formulation for chance-constrained optimization of SDLCS.

  2. 2.

We compare our proposed approach with several previously proposed techniques
and demonstrate that our method outperforms the recent techniques in [13, 14].

The proposed algorithm is demonstrated on several manipulation systems with
linear dynamics.

##  II Related Work

In this section, we review some of the work which is most close to the work
presented in this paper. Our work is closely related to TO techniques for
contact-rich systems. Contact-implicit TO techniques are becoming very popular
for performing TO for contact-rich systems, and several techniques have been
proposed for manipulation as well as legged locomotion [6, 8, 16]. All the
above techniques assume perfect model knowledge and do not consider
uncertainty.

There has been some recent work on robust TO in SNCS [13, 14, 17]. In [13],
the authors have utilized the formulation of ERM for robust TO. ERM, first
introduced in [12] for Stochastic Linear Complementarity Problem (SLCP), aims
at minimizing the expected error in satisfying the SLCP. In [13], authors use
ERM as an additional penalty term in their TO problem. However, such a
formulation does not consider the stochastic state evolution of the system
during optimization. A chance-constrained formulation for SNCS is presented in
[14]. This method augments the ERM-augmented objective in [13] with additional
chance constraints on satisfying the complementarity constraints. The
formulation ignores the stochastic evolution of system state during
optimization, and thus borrows the limitations of [13]. Furthermore, this
formulation is incapable of enforcing a constraint violation probability
smaller than $0.5$ for any degree of uncertainty. Consequently, this method is
very fragile for trajectories with the horizon, $N>1$ as the chance of
violating the constraints for such trajectories is $0.5N\geq 1$ [18]. Our
formulation addresses these weaknesses under certain simplifying assumptions
for SDLCS.

Another line of work which is relevant to understand some of our proposed work
is related to chance-constrained optimization (CCO). This has been extensively
studied in robotics as well as optimization literature[19, 20, 21]. In [19],
authors have proposed stochastic optimization formulation for open-loop
collision avoidance problems using chance constraints under Gaussian noise.
[21] uses statistical moments of the distribution to handle non-Gaussian
chance constraints. An important point to note here is that in all CC
formulation for dynamic optimization, one needs to consider the CDF function
for the joint probability distribution of all variables. However, such
distribution is generally extremely challenging to compute. Thus, in general,
the joint chance constraint is decomposed into individual chance constraints
using Boole’s inequality (see [19, 21]), which results in very conservative
approximation of the individual constraints. Our formulation utilizes Boole’s
inequality to convert the original computationally intractable joint chance
constraints into conservative but tractable independent chance constraints.

##  III Problem Preliminary

For the completeness of the paper, we provide a brief introduction to linear
complementarity problem and its stochastic form. This is followed by a problem
formulation for robust trajectory optimization for linear dynamical systems
with stochastic complementarity solution. We also point several key
differences of our approach from previous attempts for robust trajectory
optimization for stochastic complementarity system.

###  III-A Discrete-time Linear Complementarity System (DLCS)

A DLCS is a discrete-time linear dynamical system with complementarity
constraints [11] given by:

$\displaystyle x_{k+1}=Ax_{k}+Bu_{k}+C\lambda_{k+1}+g_{k}$ (1a) $\displaystyle
0\leq\lambda_{k+1}\perp Dx_{k}+Eu_{k}+F\lambda_{k+1}+h_{k}\geq 0$ (1b)

where $k$ is the time-step index, $x_{k}\in\mathbb{R}^{n_{x}}$ is the state,
$u_{k}\in\mathbb{R}^{n_{u}}$ is the control input, and
$\lambda_{k}\in\mathbb{R}^{n_{c}}$ is the algebraic variable (e.g., contact
forces). In addition, $A\in\mathbb{R}^{n_{x}\times n_{x}}$,
$B\in\mathbb{R}^{n_{x}\times n_{u}}$, $C\in\mathbb{R}^{n_{x}\times n_{c}}$,
$g_{k}\in\mathbb{R}^{n_{x}}$, $D\in\mathbb{R}^{n_{c}\times n_{x}}$,
$E\in\mathbb{R}^{n_{c}\times n_{u}}$, $F\in\mathbb{R}^{n_{c}\times n_{c}}$,
and $h_{k}\in\mathbb{R}^{n_{c}}$. The $i$-th element of vector $p_{k}$
($p_{k}$ can be $x_{k},u_{k},\lambda_{k}$) is represented as $p_{k,i}$. The
$i$-th diagonal element of matrix $P_{k}$ is represented as $P_{k,ii}$. The
notation $0\leq a\perp b\geq 0$ denotes the complementarity constraints $a\geq
0,b\geq 0,ab=0$.

Given a $x_{k},u_{k}$, an unique solution $\lambda_{k+1}$ to (1b) exists if
$F$ is P-matrix [22]. If $F$ does not satisfy the P-matrix property, it is
possible that $\lambda_{k+1}$ satisfying (1b) is non-unique or non-existent.

###  III-B Contact-Implicit Trajectory Optimization

A contact-implicit trajectory optimization for the DLCS can be formulated as:

$\displaystyle\min_{x,u,\lambda}$
$\displaystyle\;\sum_{k=0}^{N-1}J(x_{k},u_{k},\lambda_{k})$ (2a) s. t.
$\displaystyle\;x_{k+1}=Ax_{k}+Bu_{k}+C\lambda_{k+1}+g_{k},$ (2b)
$\displaystyle\;0\leq\lambda_{k+1}\perp Dx_{k}+Eu_{k}+F\lambda_{k+1}+h_{k}\geq
0,$ (2c) $\displaystyle\;\quad
x_{0}=x_{s},x_{N}=x_{g},x_{k}\in\mathcal{X},u_{k}\in\mathcal{U},\lambda_{k}\leq\lambda_{u}$
(2d)

where $x_{s},x_{g}$ represent the initial and the terminal values,
respectively, $\mathcal{X}\subseteq\mathbb{R}^{n_{x}}$ and
$\mathcal{U}\subseteq\mathbb{R}^{n_{u}}$ are convex polytopes consisting of a
finite number of linear inequality constraints, $\lambda_{u}$ is the upper
bound of $\lambda_{k}$, and $N$ is the time horizon.

While (2) is widely used in various robotic applications, it can be fragile
under uncertainty, which is often the case in model-based manipulation. Hence,
we consider a novel formulation of (2) so that the generated trajectory from
the optimization would be robust under uncertainty.

###  III-C Stochastic Discrete-time Linear Complementarity Systems (SDLCS)

We consider the following SDLCS, i.e. DLCS with uncertainty:

$\displaystyle x_{k+1}=Ax_{k}+Bu_{k}+C\lambda_{k+1}+g_{k}+w_{k}$ (3a)
$\displaystyle 0\leq\lambda_{k+1}\perp y_{k+1}\geq 0$ (3b)

where $y_{k+1}=Dx_{k}+Eu_{k}+F\lambda_{k+1}+h_{k}+v_{k}$.
$w_{k}\in\mathbb{R}^{n_{x}},v_{k}\in\mathbb{R}^{n_{c}}$ are known additive
uncertainty. We consider the case where the coefficient matrix $C$ in (3a) and
$F$ in (3b) are stochastic matrices to discuss a more realistic stochastic
effect due to complementarity constraints. This corresponds to the case when
one might have uncertainty arising from parameter identification leading to a
SDLCS. An alternative to this is to allow the complementarity variable
$\lambda_{k+1}$ to be stochastic. However, such treatment is out of the scope
of the current work. Our treatment of SDLCS leads to stochastic evolution of
system states $x_{k}$, while we treat $\lambda_{k+1}$ as deterministic. The
assumption of determinacy in $\lambda_{k+1}$ is similar to several previous
works [12], [13], [14], [23].

The authors in [13] use ERM to solve TO of SDLCS and use the following cost
function:

$\sum_{k=0}^{N-1}\left(J\left(x_{k},u_{k},\lambda_{k+1}\right)+\beta\mathbb{E}\left[\left\|\psi\left(\lambda_{k+1},y_{k+1}\right)\right\|^{2}\right]\right)$
(4)

where $\psi$ is an Nonlinear Complementarity Problem (NCP) function, $\beta$
is a weighting scalar. We compare the robustness of our formulation with (4)
in Sec V.

##  IV Robust Trajectory Optimization for SDLCS

In this section, we describe our formulation for robust TO of SDLCS. We
consider:

$\displaystyle\min_{x,u,\lambda}$
$\displaystyle\;\mathbb{E}\left[\sum_{k=0}^{N-1}J(x_{k},u_{k},\lambda_{k})\right]$
(5a) s. t. $\displaystyle\;x_{k+1}=Ax_{k}+Bu_{k}+C\lambda_{k+1}+g_{k}+w_{k},$
(5b) $\displaystyle\;\text{Pr}\left(0\leq\lambda_{k+1}\perp y_{k+1}\geq
0,x_{k}\in\mathcal{X},\forall k\right)\geq 1-\Delta,$ (5c)
$\displaystyle\;x_{0}\sim\mathcal{N}\left(x_{s},\Sigma_{s}\right),u_{k}\in\mathcal{U},\lambda_{k}\leq\lambda_{u}$
(5d)

where Pr denotes the probability of an event and $\Delta\in\left(0,0.5\right]$
is the user-defined maximum violation probability, where the probability of
violating constraints is bounded by $\Delta$. $x_{s},\Sigma_{s}$ are the mean
and covariance matrix of the state at $k=0$. $\mathcal{X}$ and $\mathcal{U}$
are convex polytopes, consisting of a finite number of linear inequality
constraints. In Sec IV-A, we describe how we convert (5c) to a tractable
optimization problem.

For clarity of presentation, we explain the reasoning behind our formulation
shown in (5). Since the underlying SDLCS is uncertain (also see Fig. 1), we
consider a chance-constrained formulation for optimization  to capture
stochastic evolution of states (see discussion in Sec I) where we impose
multiple constraints simultaneously. This is represented as joint chance
constraints for the complementarity constraints as well as the states, which
is succinctly written in Equation (5c)). Note that we represent the chance
constraints on all the variables jointly (as is common in stochastic
optimization for dynamic systems) using the cdf for the state as well as
complementarity variables. We show in the rest of this section how the joint
constraints can be decomposed into individual chance constraints using Boole’s
inequality. It is also important to note that unlike (5), the method in [13,
14] fails to capture the stochastic evolution of states in their formulation.

In this work, we make the following assumptions for (5):

  1. 1.

Noise terms $w_{k}$, $v_{k}$ follow Gaussian distribution.

  2. 2.

The complementarity variable $\lambda_{k+1}$ is deterministic.

  3. 3.

Each element of vectors $C\lambda_{k+1}$ and $F\lambda_{k+1}$ are independent
Gaussian variables.

We explain the rationale for above assumptions in Sec IV-B.

###  IV-A Joint Linear Chance Constraints

We consider joint chance constraint such that multiple constraints are
satisfied simultaneously with a prespecified probability. More specifically,
we consider the joint chance constraint (5c) so that the complementarity
constraints and state bound constraints over the whole time horizon of the
optimized trajectory are satisfied with probability $1-\Delta$. We denote the
complementarity relationship in (3b) succinctly as
$(\lambda_{k+1,i},y_{k+1,i})\in\mathcal{S}$ for $i=1,\ldots,n_{c}$. Hence, in
this optimization problem (5), we have the following joint chance constraints:

$\begin{split}&\text{Pr}\left(0\leq\lambda_{k+1}\perp y_{k+1}\geq
0,x_{k}\in\mathcal{X},\forall k\right)\geq 1-\Delta\Longleftrightarrow\\\
&\text{Pr}\left(\bigwedge_{k=0}^{N}\left(\bigwedge_{i=1}^{n_{c}}(\lambda_{k+1,i},y_{k+1,i})\in\mathcal{S}\right)\bigwedge\left(\bigwedge_{l=1}^{L}a_{l}^{\top}x_{k}\leq
b_{l}\right)\right)\\\ &\geq 1-\Delta\end{split}$ (6)

where $\bigwedge$ is the logical AND operator. $L$ represents the number of
chance constraints involving $x$ at $k$, except for the complementarity
constraints. $a_{l}\in\mathbb{R}^{n_{x}}$ is the constant vector and $b_{l}$
is a scalar.

Obtaining a cumulative distribution function (cdf) of (6) is challenging
because the joint probability of states and complementarity variables is
considered. The only way to decompose joint chance constraints is Boole’s
inequality [18] that converts the original computationally intractable joint
chance constraints into conservative but tractable independent constraints.
Hence, similar to previous works, we employ Boole’s inequality [18] to get the
conservative approximation of (6) as follows:

$\begin{split}\text{Pr}\left(\bigwedge_{k=0}^{N}\left(\bigwedge_{i=1}^{n_{c}}(\lambda_{k+1,i},y_{k+1,i})\in\mathcal{S}\right)\right)\geq
1-\Delta_{1},\\\
\text{Pr}\left(\bigwedge_{k=0}^{N}\left(\bigwedge_{l=1}^{L}a_{l}^{\top}x_{k}\leq
b_{l}\right)\right)\geq
1-\Delta_{2},\Delta_{1}=\Delta_{2}=\frac{\Delta}{2}\end{split}$ (7)

Using Boole’s inequality again, we can further obtain the conservative chance
constraints given by:

$\displaystyle\text{Pr}\left((\lambda_{k+1,i},y_{k+1,i})\in\mathcal{S}\right)\geq
1-\frac{\Delta_{1}}{Nn_{c}},$ (8a)
$\displaystyle\text{Pr}\left(a_{l}^{\top}x_{k}\leq b_{l}\right)\geq
1-\frac{\Delta_{2}}{NL},\Delta_{1}=\Delta_{2}=\frac{\Delta}{2}$ (8b)

We discuss how to handle (8a) in Sec IV-B. We formulate (8b) as its equivalent
deterministic form:

$\displaystyle\text{Pr}\left(a_{l}^{\top}x_{k}\leq b_{l}\right)\geq
1-\frac{\Delta_{2}}{NL}\Longleftrightarrow$ (9a) $\displaystyle
a_{l}^{\top}\bar{x}_{k}\leq
b_{l}-\sqrt{a_{l}^{\top}\Sigma_{x_{k}}a_{l}}{\Phi}^{-1}(1-\frac{\Delta_{2}}{NL})$
(9b)

where $\bar{x}_{k},\Sigma_{x_{k}}$ are the mean and covariance matrix of
$x_{k}$, respectively. $\Phi^{-1}$ is an inverse of the cdf of the standard
normal distribution.

###  IV-B Chance Complementarity Constraints (CCC) for SDLCS

We make the assumptions as specified in Sec IV. While a more general
formulation could allow the complementarity variable $\lambda_{k+1}$ to be
stochastic, we do not consider it here. However, we believe that allowing $C$
and $F$ to be stochastic can achieve a similar effect in SDLCS. Furthermore,
in cases where the distribution of $\lambda_{k+1}$ is known, our proposed
formulation can be easily extended to incorporate stochasticity in
$\lambda_{k+1}$. However, for brevity, we skip these details. The Gaussian
assumption on uncertainty is made primarily to allow equivalent reformulation
of the chance constraints to deterministic inequalities.

While [14] proposed a promising CCC, their formulation possesses empty
solutions when $\Delta\leq 0.5$ (see [14]). This can result in a very fragile
trajectory since the total violation probability over $N$ steps would be
always more than 1 if $N\geq 1$ (using Boole’s inequality). This is because
they use a Non-Linear Programming (NLP) formulation which needs to impose all
CCC constraints simultaneously which compete with each other.

Figure 2: Deterministic and stochastic complementarity constraints. We have
the complementarity constraints $0\leq\lambda_{k+1,i}\perp y_{k+1,i}\geq 0$
where $y_{k+1,i}$ has uncertainty and accepts the violation of $\epsilon$.

In our formulation, we decompose stochastic complementarity constraints into
two modes (see Fig. 2) as follows:

$\displaystyle\text{Pr}\left((\lambda_{k+1,i},y_{k+1,i})\in\mathcal{S}\right)\geq
1-\theta$ (10a)
$\displaystyle\Longleftrightarrow\text{Pr}\left(\begin{aligned}
\left(\lambda_{k+1,i}\geq 0,y_{k+1,i}=0\right)\\\
\bigvee\left(\lambda_{k+1,i}=0,y_{k+1,i}\geq 0\right)\end{aligned}\right)\geq
1-\theta$ (10b) $\displaystyle\Longleftrightarrow\left\\{\begin{aligned}
\lambda_{k+1,i}\geq 0,\text{Pr}\left(y_{k+1,i}=0\right)\geq 1-\theta\\\
\text{or }\lambda_{k+1,i}=0,\text{Pr}\left(y_{k+1,i}\geq 0\right)\geq
1-\theta\end{aligned}\right.$ (10c)

where $\theta=\frac{\Delta_{1}}{Nn_{c}}$. Note that now
$y_{k+1}\sim\mathcal{N}\left(\bar{y}_{k+1},\Sigma_{y_{k+1}}\right)$. To
realize lower violation probabilities, we first propose the following CCC
using MIP as:

$\displaystyle z_{k,i,0}\geq 0,\Longrightarrow\lambda_{k+1,i}\geq
0,\text{Pr}\left(y_{k+1,i}=0\right)\geq 1-\theta,$ (11a) $\displaystyle
z_{k,i,1}\geq 0,\Longrightarrow\lambda_{k+1,i}=0,\text{Pr}\left(y_{k+1,i}\geq
0\right)\geq 1-\theta$ (11b)

where $z_{k,i,0},z_{k,i,1}$ represent the integer variables to represent the
two modes which satisfies $z_{k,i,0}+z_{k,i,1}=1$ for $i$-th complementarity
constraint at instant $k$.

However, $\text{Pr}\left(y_{k+1,i}=0\right)$ is zero (as probability measure
for singleton sets is zero) so that we cannot directly use (11). To alleviate
this issue while avoiding negative values for $\lambda$, we propose the
following CCC using a relaxation for complementarity constraints (see Fig. 2):

$\displaystyle z_{k,i,0}\geq 0,\Longrightarrow\lambda_{k+1,i}\geq
0,\text{Pr}\left(0\leq y_{k+1,i}\leq\epsilon\right)\geq 1-\theta,$ (12a)
$\displaystyle z_{k,i,1}\geq
0,\Longrightarrow\lambda_{k+1,i}=0,\text{Pr}\left(y_{k+1,i}\geq\epsilon\right)\geq
1-\theta$ (12b)

where $\epsilon>0$ is the acceptable violation in the complementarity
constraints.

We have two-sided linear chance constraints in (12a). We decompose (12a) as
two one-sided chance constraints so that we can use the same reformulation in
(9). Note that each one-sided chance constraints, obtained from the two-sided
chance constraint, are formulated with a maximum violation probability of
$\frac{\theta}{2}$.

Since we have integer constraints, MIP can impose individual constraints for
each mode. Thus, we do not need to impose all constraints simultaneously like
the NLP formulation in [14]. This provides a lower bound for $\theta$ as
function of $\epsilon$, $\bar{y}_{k+1,i}$, and $\Sigma_{y_{k+1},ii}$, which is
presented as a lemma.

######  Lemma IV.1.

Suppose the CCC are formulated as (12) and $\epsilon$, $\bar{y}_{k+1,i}$, and
$\Sigma_{y_{k+1},ii}$ are specified. Then (i) (12a) is feasible for all
$\theta>2(1-\Phi(\frac{\epsilon}{2\Sigma_{y_{k+1},ii}}))$ and (ii) (12b) is
feasible for all
$\theta>1-\Phi\left((\bar{y}_{k+1,i}-\epsilon)/\Sigma_{y_{k+1}ii}\right)$.

###### Proof.

Consider case (i): From (9b) and (12a), the two-side chance constraints in
(12a) are converted to their deterministic
forms:$\Sigma_{y_{k+1},ii}\Phi^{-1}\left(1-\theta/2\right)\leq\bar{y}_{k+1,i}\leq\epsilon-\Sigma_{y_{k+1},ii}\Phi^{-1}\left(1-\theta/2\right)$.
To have a nonempty solution,
$\epsilon-2\Sigma_{y_{k+1},ii}\Phi^{-1}\left(1-\theta/2\right)>0$. Simplifying
this equation, we obtain the bound specified in (i). Consider case (ii): From
(9b) and (12b), the one-side chance constraints in (12b) are converted as:
$\bar{y}_{k+1,i}\geq\epsilon+\Sigma_{y_{k+1},ii}\Phi^{-1}\left(1-\theta\right)$.
Simplifying this equation, we obtain the bound specified in (ii). ∎

Remark 1: From Lemma IV.1, it is easy to show that $\theta<\frac{1}{2}$ if
$\frac{\epsilon}{2\Sigma_{y_{k+1},ii}}>\Phi^{-1}(\frac{3}{4})$ for case (i),
and if $(\bar{y}_{k+1,i}-\epsilon)/\Sigma_{y_{k+1}ii}>\Phi^{-1}(\frac{1}{2})$
for case (ii). In contrast, the formulation in [14] cannot enforce the chance
constraints for any $\theta<0.5$.

We use the following equations for uncertainty propagation in the SDLCS:

$\displaystyle\bar{x}_{k+1}=A\bar{x}_{k}+Bu_{k}+\overline{C\lambda}_{k+1}+g_{k}+\bar{w}_{k},$
(13a)
$\displaystyle\Sigma_{x_{k+1}}=A\Sigma_{x_{k}}A^{\top}+\Sigma_{C\lambda_{k+1}}+W$
(13b)

where $W$ represents the noise covariance matrix and
$\overline{C\lambda}_{k+1}$ represents a mean of ${C\lambda}_{k+1}$.
$\Sigma_{C\lambda_{k+1}}=\mathbb{E}\left[\left(C\lambda_{k+1}-\overline{C\lambda}_{k+1}\right)\left({C\lambda_{k+1}-\overline{C\lambda}_{k+1}}\right)^{\top}\right]$,
which is a diagonal matrix because of the independence of random variables.
Note that $\lambda_{k+1}$ is a decision variable. Consequently, we introduce
another simplification by considering the worst-case uncertainty for
$\lambda_{k+1}$ during uncertainty propagation. This conservative
simplification offers computational advantages during the resulting
optimization.

###  IV-C Mixed-Integer Quadratic Programming with Chance Constraints (MIQPCC)

In this section, we present our MIQPCC formulation to solve (5). To impose our
proposed CCC, one can solve either MIP or NLP. Our MIP-based method solves
disjunctive inequalities while NLP needs to impose all CCC simultaneously,
which yields an empty solution for $\Delta\leq 0.5$.

Our proposed MIQPCC is formulated as follows:

$\displaystyle\min_{x,u,\lambda,z}\sum_{k=0}^{N-1}\bar{x}_{k}^{\top}Q\bar{x}_{k}+u_{k}^{\top}Ru_{k}$
(14a) $\displaystyle\text{s. t.
}\bar{x}_{k+1}=A\bar{x}_{k}+Bu_{k}+\overline{C\lambda}_{k+1}+g_{k}+\bar{w}_{k},$
(14b)
$\displaystyle\Sigma_{x_{k+1}}=A\Sigma_{x_{k}}A^{\top}+\Sigma_{w,C\lambda_{k+1}}+W,$
(14c) $\displaystyle
x_{0}\sim\mathcal{N}\left(x_{s},\Sigma_{s}\right),u_{k}\in\mathcal{U},\lambda_{k}\leq\lambda_{u},$
(14d) $\displaystyle a_{l}^{\top}\bar{x}_{k}\leq b_{l}-\alpha\kappa,$ (14e)
$\displaystyle z_{k,i,0}+z_{k,i,1}=1,$ (14f) $\displaystyle
0\leq\lambda_{k+1,i}\leq Mz_{k,i,0},$ (14g) $\displaystyle\zeta\psi
z_{i,k,0}+(\epsilon+\eta\psi)z_{k,i,1}\leq\bar{y}_{k+1,i}$ (14h)
$\displaystyle\bar{y}_{k+1,i}\leq(\epsilon-\zeta\psi)z_{k,i,0}+Mz_{k,i,1},$
(14i)

where $Q=Q^{\top}\geq 0,R=R^{\top}>0$,
$\alpha={\Phi}^{-1}(1-\frac{\Delta}{2NL}),\zeta={\Phi}^{-1}(1-\frac{\Delta}{4Nn_{c}}),\eta={\Phi}^{-1}(1-\frac{\Delta}{2Nn_{c}}),\kappa=\sqrt{a_{l}^{\top}\Sigma_{x_{k}}a_{l}},\psi=\sqrt{\Sigma_{y_{k+1,ii}}}$.
$z_{k,i,0},z_{k,i,1}$ are the binary decision variables for the $i$-th
complementarity constraint at $k$ to represent mode 1, 2, respectively. Using
these binary variables, we employ big-M formulation to deal with disjunctive
inequalities in our CCC. The parameter $M$ is a valid upper bound for
$\lambda_{k},y_{k}$.

##  V Numerical Simulations

We validate our proposed methods for three benchmark DLCS: a cartpole with
softwalls, a sliding box with friction, and dual manipulators as illustrated
in Fig. 3, inspired by [10]. Through the experiments, we try to answer the
following questions:

  1. 1.

Can our proposed optimization generate robust open-loop trajectories?

  2. 2.

Can our proposed formulation satisfy the probabilistic constraints imposed
during optimization?

  3. 3.

How does the proposed method compare against the previous methods for robust
optimization in SDLCS?

###  V-A Implementation Details

We implemented our method in Python using Gurobi [24] to solve the proposed
MIQP. We implemented the MPCC with PYROBOCOP [6] to solve the ERM-based method
in [13] and the CCC method in [14]. The examples are implemented on a computer
with Intel i7-8565U processor.

To verify the robustness of open-loop trajectories obtained from our proposed
optimization, we use Monte Carlo simulations. We propagate the dynamics by
finding the roots of the complementarity system with sampled parameters given
the control sequence obtained from optimization. We run each case for 1000
trials with different sampled parameters to estimate the probability of
failure. Note that unlike the continuous-domain dynamics, we cannot rollout
the dynamics for SDLCS with the given control sequences since we do not have
the access to $\lambda_{k+1}$. We add the noise sampled from the distribution
which was used during optimization.

For simplicity, we show the continuous-time dynamics. We then discretize
continuous-time dynamics into discrete-time dynamics using the explicit Euler
method with sample time $dt=0.033$. For notation simplicity, we denote
$x_{0},\Sigma_{0}$ as mean and covariance matrix at $k=0$ for states of
systems, respectively.

###  V-B Example Details

####  V-B1 Cartpole with Softwalls

The continuous-time dynamics with complementarity constraints for the cartpole
with softwalls is as follows:

$\displaystyle\dot{x}_{1}=x_{3},\dot{x}_{2}=x_{4},\dot{x}_{3}=g\frac{m_{p}}{m_{c}}x_{2}+\frac{1}{m_{c}}u_{1},$
(15a)
$\displaystyle\dot{x}_{4}=\frac{g\left(m_{c}+m_{p}\right)}{lm_{c}}x_{2}+\frac{1}{lm_{c}}u_{1}+\frac{\lambda_{1}}{lm_{p}}-\frac{\lambda_{2}}{lm_{p}},$
(15b) $\displaystyle 0\leq\lambda_{1}\perp
lx_{2}-x_{1}+\frac{1}{k_{1}}\lambda_{1}+d\geq 0,$ (15c) $\displaystyle
0\leq\lambda_{2}\perp x_{1}-lx_{2}+\frac{1}{k_{2}}\lambda_{2}+d\geq 0$ (15d)

where $x_{1}$ is the cart position, $x_{2}$ is the pole angle, the $x_{3}$ and
$x_{4}$ are their derivatives. $u_{1}$ is the control and
$\lambda_{1},\lambda_{2}$ are the reaction forces at from the wall 1, 2,
respectively. We consider the additive noise $w$, the zero-mean i.i.d.
Gaussian noise which standard deviation is $2\times 10^{-4}$, to
$x_{1,k},x_{2,k}$. $k_{1}=10,k_{2}=10$ are the stiffness of walls $1$ and $2$,
respectively. In this example, we assume that the uncertainty also arises from
the $\frac{1}{k_{1}},\frac{1}{k_{2}}$ which standard deviations are $10^{-5}$.
$g=9.81$ is the gravitational acceleration, $m_{p}=0.1,m_{c}=1.0$ are the mass
of the pole, cart, respectively. $l=0.5$ is the length of the pole and
$d=0.15$ is the distance from the origin of the coordinate to the walls.

The optimization setup is as follows.
$N=20,M=100,Q=\text{diag}(0,0,0,0),R=0.01$, $\epsilon=0.002,$
$x_{0}=[-0.15,0,0,0]^{\top},\Sigma_{0}=\text{diag}(0,0,0,0)$. We also have the
following chance constraints: $\text{Pr}(x_{1,k}\leq 0.05)\geq
1-\frac{\Delta}{4N},\text{Pr}(x_{2,k}\leq 0.15)\geq
1-\frac{\Delta}{4N},\forall k=0,\ldots,N-1$, $\text{Pr}(-0.02\leq x_{1,N}\leq
0.02)\geq 1-\frac{\Delta}{4N},\text{Pr}(-0.04\leq x_{2,N}\leq 0.04)\geq
1-\frac{\Delta}{4N}$.

Figure 3: Problems described in Sec V. (a) cartpole with softwalls, (b)sliding
box with friction, and (c)dual manipulators.

####  V-B2 Sliding Box with Friction

The continuous-time quasi-static dynamics with complementarity constraints for
sliding box with Coulomb friction is as follows:

$\displaystyle\dot{x}_{1}=x_{2},\alpha\dot{{x}}_{1}=u+\lambda_{+}-\lambda_{-},$
(16a) $\displaystyle 0\leq\gamma\perp\mu mg-\lambda_{+}-\lambda_{-}\geq 0,$
(16b) $\displaystyle 0\leq\lambda_{+}\perp\gamma+u+\lambda_{+}-\lambda_{-}\geq
0,$ (16c) $\displaystyle
0\leq\lambda_{-}\perp\gamma-u-\lambda_{+}+\lambda_{-}\geq 0$ (16d)

$x_{1}$ is the box position and $x_{2}$ is the box velocity. $u$ is the
control and $\lambda_{+},\lambda_{-}$ are the positive and negative component
of the friction force, respectively. $\gamma$ is the slack variable.
$\alpha=4$ is the damping constant, $m=1$ is the mass of the box, and
$\mu=0.1$ is the coefficient of friction. We consider the additive i.i.d.
Gaussian noise $w$ as ${x_{1,k+1}}=x_{1,k}+x_{2,k}dt+w$. The standard
deviations of $w$ is $4\times 10^{-4}$. $g=9.81$ is the gravitational
acceleration. We assume that the uncertainty also arises from the $\mu$ which
standard deviations are $10^{-5}$.

The optimization setup is as follows.
$N=20,M=100,Q=\text{diag}(0,0,0,0),R=0.01$,
$\epsilon=0.01,x_{0}=[1,-1]^{\top},\Sigma_{0}=\text{diag}(0,0)$. We also have
the following chance constraints: $\text{Pr}(x_{1,k}\geq 0.885)\geq
1-\frac{\Delta}{2N},\forall k=0,\ldots,N-1$, $\text{Pr}(0.89\leq x_{1,N}\leq
0.91)\geq 1-\frac{\Delta}{2N},\text{Pr}(-0.1\leq x_{2,N}\leq 0.1)\geq
1-\frac{\Delta}{2N}$.

####  V-B3 Dual Manipulators

We consider the example where the box is manipulated by two manipulators with
Coulomb friction and the contact forces from the manipulators. The continuous-
time quasi-static dynamics is as follows:

$\displaystyle\dot{x}_{1}=x_{2},\alpha\dot{x}_{1}=\lambda_{1}-\lambda_{2}+\lambda_{+}-\lambda_{-},$
(17)
$\displaystyle\dot{x}_{3}=x_{4},\dot{x}_{4}=u_{1},\dot{x}_{5}=x_{6},\dot{x}_{6}=u_{2},$
$\displaystyle 0\leq\lambda_{1}\perp x_{1}-x_{3}+\frac{1}{k}\lambda_{1}\geq
0,$ $\displaystyle 0\leq\lambda_{2}\perp
x_{5}-x_{1}+\frac{1}{k}\lambda_{2}\geq 0,$ $\displaystyle 0\leq\gamma\perp\mu
mg-\lambda_{+}-\lambda_{-}\geq 0,$ $\displaystyle
0\leq\lambda_{+}\perp\gamma+\lambda_{1}-\lambda_{2}+\lambda_{+}-\lambda_{-}\geq
0,$ $\displaystyle
0\leq\lambda_{-}\perp\gamma-\lambda_{1}+\lambda_{2}-\lambda_{+}+\lambda_{-}\geq
0$

$x_{1},x_{3},x_{5}$ are the positions of the box, the left arm, the right arm,
respectively and $x_{2},x_{4},x_{6}$ are their derivatives. $u_{1},u_{2}$
represent the controls of the left and the right arm, respectively.
$\lambda_{+},\lambda_{-}$ are the positive and negative component of the
friction force, respectively. $\gamma$ is the slack variable.
$\lambda_{1},\lambda_{2}$ are the contact forces from the left arm and the
right arm, respectively. We set $g=9.81$, $m=1$, $k=100$, $\mu=0.1$. We
discretize the dynamics (17) with $dt=0.033$ and add the zero-mean i.i.d.
Gaussian noise $w$ which standard deviation is 0.0002 such as
${x_{1,k+1}}=x_{1,k}+x_{2,k}dt+w$. The standard derivation of $\mu$ and
$\frac{1}{k}$ are 0.0001.

The optimization setup in this example is as follows.
$N=20,M=50,Q=\text{diag}(0,0,0,0,0,0),R=\text{diag}(1,1),\epsilon=0.0042$,
$x_{0}=[0.1,-1.1,0,0,0.1,0]^{\top},\Sigma_{0}=\text{diag}(0,0,0,0,0,0)$. We
have the following chance constraints: $\text{Pr}(x_{1,k}\geq-0.17)\geq
1-\frac{\Delta}{2N},\forall k=0,\ldots,N-1$, $\text{Pr}(-0.01\leq x_{1,N}\leq
0.01)\geq 1-\frac{\Delta}{2N}$.

###  V-C Robustness of Open-Loop Trajectories

The optimized control and state trajectories for the three systems using our
proposed method are shown in Fig. 4-Fig. 6. Overall, these figures show that
the planner generates state trajectories that are farther away from the bound
specified in the chance constraints as the violation probability decreases.
For instance, Fig. 4 shows that the trajectories are farther away from
$x=0.05$ as $\Delta$ decreases. In addition, the trajectory with $\Delta=0.02$
reaches its maximum value earlier than other trajectories to account for the
evolution of the uncertainty. We observe the same behavior for the other
example too. In addition, these figures illustrate that the control costs
increase as $\Delta$ decreases. This illustrates the trade-off relation
between safety and cost.

At this point, we would like to discuss the magnitude of uncertainty we
consider in these problems. Compared to some other stochastic optimal control
works [19, 20], the uncertainty in these problems is relatively smaller. There
are several reasons why we need to have a smaller uncertainty. Note that as we
have explained in Sec III, our approach satisfies joint constraints on
multiple constraints together. First, our formulation has chance
complementarity constraints in addition to chance constraints on states, which
are commonly used. Our formulation has more number of chance constraints, and
consequently, the lower uncertainty is required because of the conservative
approximation of Boole’s inequality to resolve joint chance constraints into
individual constraints as explained in Sec IV-A, Sec IV-B. Second, we need to
have a small $\epsilon$ to avoid large violation of complementarity
constraints, which requires small uncertainty. Finally, we would like to
emphasize that allowing larger uncertainties requires either better resolution
of joint chance constraints or covariance steering approaches [20], which is
out of scope for the current study.

Figure 4: State and control trajectories with different $\Delta$ for the
cartpole example. First, the cart moves in the negative direction to utilize
the contact force $\lambda_{2}$ because the control input is bounded. Once the
cart obtains enough $\lambda_{2}$, the cart is accelerated to the positive
direction. We can observe the effect of our proposed chance constraints in
particular around $t\in[0,0.1]$ and $t\in[0.4,0.5]$. When $t\in[0,0.1]$, the
mode changes from the ”contact on the wall 2” to the ”no contact” and the cart
tries to be far from wall 2 to satisfy the CCC. When $t\in[0.4,0.5]$, the
trajectories are farther away from $x_{1}=0.05$ and $x_{2}=0.15$ as $\Delta$
decreases. Figure 5: State and control trajectories with different $\Delta$
for a sliding box with friction. First, the box is accelerated in the positive
direction. Then, the control decreases with time to regulate the box around
the origin by employing the friction forces. We can observe the effect of our
proposed chance constraints in particular around $t\in[0.2,0.3]$ where the
trajectories are farther away from $x_{1}=0.88$ as $\Delta$ decreases. Figure
6: Time history of $x_{1}$ with different maximum violation probabilities
$\Delta$ for dual manipulation. First, the box is pushed by the right arm in
the negative direction. Next, the left arm regulates the box to the origin. In
particular, around $t\in 0.2-0.3$ s, the trajectories are farther away from
$x_{1}=-0.17$ as $\Delta$ decreases.

###  V-D Monte Carlo Simulation Results

Table I-Table III show our Monte Carlo simulation results given the control
sequences with different $\Delta$ from the optimization compared to the ERM
method in [13] and the CCC method in [14]. We run the ERM method in [13] with
different weighting $\beta$ and the CCC [14] with violation probability
$\Delta_{z}=0.5$. $\beta$ was chosen so that the magnitude of the ERM cost is
a similar order of other costs. For a fair comparison, we regard that the
constraints are violated if the chance constraints are not satisfied in our
method. We regard that the constraints are violated in the ERM in [13] and the
CCC methods in [14] if the terminal chance constraints used in our proposed
method are not satisfied.

Table I: Comparison of our specified $\Delta$ in optimization, specified
$\beta$ in ERM in [13], and the CCC in [14] with $\Delta_{z}=0.5$, and
obtained $\Delta$ from the simulation of ”cartpole with softwalls” over 1000
samples.

| $\Delta=0.5$ | $\Delta=0.2$ | $\Delta=0.1$ | $\Delta=0.02$
---|---|---|---|---
Obtained $\Delta$ | 0.190 | 0.147 | 0.085 | 0.020

| $\beta=10^{3}$ | $\beta=10^{4}$ | $\beta=10^{5}$ | $\Delta_{z}=0.5$
---|---|---|---|---
Obtained $\Delta$ | 0.75 | 1.0 | 1.0 | 0.91

Table II: Comparison of our specified $\Delta$ in optimization, specified
$\beta$ in ERM in [13], and the CCC in [14] with $\Delta_{z}=0.5$, and
obtained $\Delta$ from the simulation of ”a sliding box with friction” over
1000 samples.

| $\Delta=0.5$ | $\Delta=0.1$ | $\Delta=0.01$ | $\Delta=0.002$
---|---|---|---|---
Obtained $\Delta$ | 0.080 | 0.051 | 0.027 | 0.010

| $\beta=10^{3}$ | $\beta=10^{4}$ | $\beta=10^{5}$ | $\Delta_{z}=0.5$
---|---|---|---|---
Obtained $\Delta$ | 1.0 | 1.0 | 1.0 | 0.91

Table III: Comparison of our specified $\Delta$ in optimization, specified
$\beta$ in ERM, and the CCC with $\Delta_{z}=0.5$, and obtained $\Delta$ from
the simulation of ”dual manipulation” over 1000 samples.

| $\Delta=0.5$ | $\Delta=0.4$ | $\Delta=0.3$ | $\Delta=0.2$
---|---|---|---|---
Obtained $\Delta$ | 0.419 | 0.317 | 0.257 | 0.217

| $\beta=10^{3}$ | $\beta=10^{4}$ | $\beta=10^{5}$ | $\Delta_{z}=0.5$
---|---|---|---|---
Obtained $\Delta$ | 1.0 | 1.0 | 1.0 | 1.0

Table I shows that the empirically obtained violation probabilities are lower
than the specified violation probabilities used in our proposed optimization.
In contrast, the control sequences based on the ERM method in [13] show
$100\%$ violation probabilities with $\beta=10^{4},10^{5}$, which are much
worse than the obtained violation probabilities using our proposed method.
With $\beta=10^{3}$, we got a relatively good violation probability. The CCC
in [14] could also show the relatively good violation probability compared to
the ERM-based method with $\beta=10^{4},10^{5}$ but shows the worse violation
probability compared to our method with $\Delta=0.5$ and the ERM with
$\beta=10^{3}$. Thus, we confirm that our proposed approach satisfies chance
constraints in the simulator in this example.

Table II shows that empirically obtained violation probabilities are lower
than the specified violation probabilities used in our proposed optimization
like the cartpole example, except for the cases with $\Delta=0.01,0.002$.
There are several factors that contribute to the violation of the chance
constraints. Unlike the cartpole example, $F$ is not a P matrix so we can get
the multiple solutions in $\lambda$, which can lead to non-Gaussian
distributions. Also, even though $\epsilon$ is small, it is not zero so the
actual trajectory in the simulator cannot be exactly the same as the
trajectory from the optimization even with no noise. While we can ignore these
effects with relatively large $\Delta$, we cannot ignore these effects anymore
with the small $\Delta$. Although the planner could not satisfy the chance
constraints for all $\Delta$ in this example, our method achieves much lower
violation probabilities compared to the ERM in [13] and the CCC in [14]. Table
III shows that we have the same discussion for the dual manipulators example
as for the pushing a box example.

Fig. 7 and Fig. 8 show that our proposed planner could successfully drive the
system to the goal state. We also observe that with decreasing $\Delta$, the
system trajectories move further away from state set boundaries to satisfy
tighter chance constraints. For Fig. 9, while the majority of the sampled
trajectories converge to the specified terminal constraints, some of them
clearly converged to other states. This result also shows that the true
distribution of the uncertainty for the dynamics systems with LCS is not
Gaussian.

Figure 7: Simulated trajectories of $x_{2}$ of the cartpole example over 1000
samples with $\Delta=0.5$ for the left column and with $\Delta=0.02$ for the
right column. The bottom row enlarges the top row figures around the area
where the chance constraints effect is observed. The red line shows the 99.9
$\%$ confidence interval. Figure 8: Simulated trajectories of $x$ of sliding
box with friction example over 1000 samples with $\Delta=0.5$ for the left
column and with $\Delta=0.002$ for the right column. The bottom row enlarges
the top row figures around the area where the chance constraints effect is
observed. The red line shows the 99 $\%$ confidence interval. Figure 9:
Simulated trajectories of $x_{1}$ over 1000 samples with $\Delta=0.2$ for dual
manipulation. The right figure shows the enlarged figure of the left figure.
The red line shows the 99 $\%$ confidence interval.

##  VI Discussion and Conclusion

The hybrid dynamics of frictional interaction as well as uncertainty
associated with frictional parameters make the efficient design of model-based
controllers for manipulation challenging. In this paper, we presented a robust
TO technique for contact-rich systems. We presented a formulation for chance
constrained optimization for SDLCS which is solved using MIQPCC. We compared
our proposed approach against other recent techniques for robust optimization
for stochastic complementarity systems. We showed that our formulation leads
to more robust trajectories compared to these techniques.

In the future, we would like to relax certain assumptions in this work. We
would like to propose solutions for general non-linear stochastic
complementarity systems in the presence of non-Gaussian noise. In the current
work, using joint chance constraints on all the variables results in
conservative solutions. To consider these problems, the study of nonlinear
uncertainty propagation in SNCS is required. Also, we need to solve mixed-
integer non-linear programming. We would also like investigate how we can
relax the conservative solutions obtained by our proposed approach using
better measures for risk. We would also like to incorporate real-time sensor
input [25] to develop algorithms for stochastic model predictive control of
complex manipulation problems [17]. Another interesting line of work would to
be to include a Reinforcement learning algorithm to get model updates [26]
during learning.

## References

  * [1] E. Todorov, “Implicit nonlinear complementarity: A new approach to contact dynamics,” in _Proc. 2010 IEEE Int. Conf. Robot. Automat._ IEEE, 2010, pp. 2322–2329.
  * [2] E. Drumwright and D. A. Shell, “An evaluation of methods for modeling contact in multibody simulation,” in _Proc. 2011 IEEE Int. Conf. Robot. Automat._ IEEE, 2011, pp. 1695–1701.
  * [3] E. Drumwright and D. Shell, “Modeling contact friction and joint friction in dynamic robotic simulation using the principle of maximum dissipation,” in _Algorithmic foundations of robotics IX_. Springer, 2010, pp. 249–266.
  * [4] M. Anitescu and F. A. Potra, “Formulating dynamic multi-rigid-body contact problems with friction as solvable linear complementarity problems,” _Nonlinear Dynamics_ , vol. 14, no. 3, pp. 231–247, 1997.
  * [5] Y. Shirai, X. Lin, Y. Tanaka, A. Mehta, and D. Hong, “Risk-aware motion planning for a limbed robot with stochastic gripping forces using nonlinear programming,” _IEEE Robot. Autom. Lett._ , vol. 5, no. 4, pp. 4994–5001, 2020.
  * [6] A. U. Raghunathan, D. K. Jha, and D. Romeres, “PYROBOCOP : Python-based robotic control & optimization package for manipulation and collision avoidance,” _CoRR_ , vol. abs/2106.03220, 2021.
  * [7] S. Jin, D. Romeres, A. Ragunathan, D. K. Jha, and M. Tomizuka, “Trajectory optimization for manipulation of deformable objects: Assembly of belt drive units,” _arXiv preprint arXiv:2106.00898_ , 2021.
  * [8] M. Posa, C. Cantu, and R. Tedrake, “A direct method for trajectory optimization of rigid bodies through contact,” _Int. J. Rob. Res._ , vol. 33, no. 1, pp. 69–81, 2014.
  * [9] M. K. Camlibel, J.-S. Pang, and J. Shen, “Lyapunov stability of complementarity and extended systems,” _SIAM J Optimization_ , vol. 17, no. 4, pp. 1056–1101, 2006.
  * [10] A. Aydinoglu, V. M. Preciado, and M. Posa, “Stabilization of complementarity systems via contact-aware controllers,” _CoRR_ , vol. abs/2008.02104, 2020\.
  * [11] A. U. Raghunathan and J. L. Linderoth, “Stability analysis of discrete-time linear complementarity systems,” _arXiv_ , 2020. [Online]. Available: arXiv:2012.13287
  * [12] X. Chen and M. Fukushima, “Expected residual minimization method for stochastic linear complementarity problems,” _Mathematics of Operations Research_ , vol. 30, no. 4, pp. 1022–1038, 2005.
  * [13] L. Drnach and Y. Zhao, “Robust trajectory optimization over uncertain terrain with stochastic complementarity,” _IEEE Robot. Autom. Lett._ , vol. 6, no. 2, pp. 1168–1175, 2021.
  * [14] L. Drnach, J. Z. Zhang, and Y. Zhao, “Mediating between contact feasibility and robustness of trajectory optimization through chance complementarity constraints,” _Front. Robo. AI_ , vol. 8, 2021.
  * [15] A. Ben-Tal, L. El Ghaoui, and A. Nemirovski, _Robust optimization_. Princeton university press, 2009.
  * [16] Z. Manchester and S. Kuindersma, “Variational contact-implicit trajectory optimization,” in _Robotics Research_. Springer, 2020, pp. 985–1000.
  * [17] Y. Shirai, D. Jha, A. Raghunathan, and D. Romeres, “Robust pivoting: Exploiting frictional stability using bilevel optimization,” in _Proc. 2022 IEEE Int. Conf. Robot. Automat._ , 2022, p. accepted.
  * [18] A. Prékopa, “Boole-bonferroni inequalities and linear programming,” _Operations Research_ , vol. 36, no. 1, pp. 145–162, 1988.
  * [19] L. Blackmore, M. Ono, and B. C. Williams, “Chance-constrained optimal path planning with obstacles,” _IEEE Transactions on Robotics_ , vol. 27, no. 6, pp. 1080–1094, 2011.
  * [20] K. Okamoto and P. Tsiotras, “Optimal stochastic vehicle path planning using covariance steering,” _IEEE Robot. Autom. Lett._ , vol. 4, no. 3, pp. 2276–2281, 2019.
  * [21] A. Wang, A. Jasour, and B. C. Williams, “Non-gaussian chance-constrained trajectory planning for autonomous vehicles under agent uncertainty,” _IEEE Robot. Autom. Lett._ , vol. 5, no. 4, pp. 6041–6048, 2020.
  * [22] R. Cottle, J. Pang, and R. Stone, _The Linear Complementarity Problem_ , ser. Classics in Applied Mathematics. Society for Industrial and Applied Mathematics, 2009.
  * [23] Y. Tassa and E. Todorov, “Stochastic complementarity for local control of discontinuous dynamics,” in _Proceedings of Robotics: Science and Systems_ , Zaragoza, Spain, June 2010.
  * [24] Gurobi Optimization, LLC, “Gurobi Optimizer Reference Manual,” 2021. [Online]. Available: https://www.gurobi.com
  * [25] S. Dong, D. Jha, D. Romeres, S. Kim, D. Nikovski, and A. Rodriguez, “Tactile-RL for insertion: Generalization to objects of unknown geometry,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , 2021.
  * [26] K. Ota, D. K. Jha, D. Romeres, J. van Baar, K. A. Smith, T. Semitsu, T. Oiki, A. Sullivan, D. Nikovski, and J. B. Tenenbaum, “Data-efficient learning for complex and real-time physical problem solving using augmented simulation,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 4241–4248, 2021\.


# A Model-based Deep Learning Reconstruction for X-ray CT

Kaichao Liang, Li Zhang, Yirong Yang, Hongkai Yang, Yuxiang Xing ∗
Key Laboratory of Particle & Radiation Imaging, Tsinghua University
Beijing, China 100084
∗Corresponding author: Yuxiang Xing Email<EMAIL_ADDRESS>

###### Abstract

Low dose CT is of great interest in these days. Dose reduction raises noise
level in projections and decrease image quality in reconstructions. Model
based image reconstruction can combine statistical noise model together with
prior knowledge into an Bayesian optimization problem so that significantly
reduce noise and artefacts. In this work, we propose a model-base deep
learning for CT reconstruction so that a reconstruction network can be trained
with no ground-truth images needed. Instead of minimizing cost function for
each image, the network learns to minimize an ensemble cost function for the
whole training set. No iteration will be needed for real data reconstruction
using such a trained network. We experimented with a penalized weighted least-
squares (PWLS) cost function for low dose CT reconstruction and tested on data
from a practical dental CT. Very encouraging results with great noise
reductions are obtained.

##  1 Introduction

As X-ray CT is a widely used medical imaging modality for clinical diagnosis
tool, its radiation dose is of big concerned. High level of dose received by
patients increase the risk of cancer [1]. In the past decades, lowering X-ray
dose by adjusting tube-voltage or tube-current-time have been studied
extensively [2, 3, 4]. However, further reducing incident photons can cause
high level statistical noise, and electronic noise becomes non-neglectable as
well. If directly reconstructed with an FBP type algorithm, the noise in
projection data will be magnified by the high-pass filtration (RL filter or SL
filter) and degrade image quality.

Previously, many researches have been conducted on low dose CT (LDCT) image
optimization[1]. These algorithms generally can be divided into three
categories. The first category is filtration method. Manually designed filters
or filtration algorithms are adopted to suspense noise. Pre-processing methods
focus on projection domain denoising [5], and post-processing methods focus on
image domain denoising after analytical reconstruction [6]. Zhang et al.
proposed a LDCT local texture restoration algorithm with priors extracted from
NDCT image dataset [7], which integrated NDCT reference dataset into
filtration methods. However, low-pass filters restrain noise but generally
cause detail information loss in projection data which further cause image
blurring. The second category is iterative reconstruction methods [8, 9, 10].
Model-based iterative reconstruction (MBIR) is based on Bayesian theory and
models noise model, physical process and system configuration into
reconstruction. For example, combining statistical iterative reconstruction
(SIR) with none-local-mean (NLM) constraint for low-dose CT reconstruction [8,
9] to suppress noise. Instead of using NLM, Xu adopted dictionary learning
which learns prior constraint based on normal-dose CT dataset [10]. The images
reconstructed from low-dose projections are expected to be of similar
distribution as normal dose CT reconstructions. However, iterative process is
very time consuming which is obstructed its usage in real-time reconstruction.
Besides, hyper-parameters in cost function and number of iteration need to be
tuned according to different data.

Recently, deep learning methods have been widely conducted in LDCT problems
[11, 12, 13, 14]. Hu and his colleagues used a three-layer CNN structure for
LDCT image optimization and achieved visible image quality improvements [11].
Later, Chen proposed a more complex encoder-decoder structure CNN for LDCT and
achieved better results [12]. Convolutional neural networks (CNN) learns the
relationship between LDCT images and NDCT images from training set. Because of
the large cost of real data, LDCT and corresponding NDCT images are often from
simulations. In these works, normal dose CT (NDCT) images are required to form
labels, which is easy to access in simulation study. However, in practical
complex situation, features of real data can be very different from
simulation, which could heavily influences the generalization of networks. No
ground-truth images available for labels in real situation become an obstacle
for practical problems.

In this work, we propose a model based deep learning (MBDL) method for CT
reconstruction that requires no ground-truth information. The network can be
directly trained with LDCT projection data only, encouraging reconstruction
results from practical low dose dental CT data are obtained.

##  2 Methods

###  2.1 Cost function to model X-ray CT imaging

In regression problem, labels are usually used to form L1/L2-norm as the cost
function of networks to minimize. In many situations, ground-truth information
might be available. The basic idea of our method is to incorporate a system
model based cost function for deep-learning network. Without losing the
generality, we took Bayesian estimation principle to construct our cost
function for CT reconstruction similar to MBIR methods. In general, Bayesian
cost function contains two terms: data fidelity term which models the noise of
projection data, and prior term which constraint the reconstructed image to
follow prior information about general image characteristics.

$\Phi\left(\bm{\upmu}\right)=\phi_{\rm
f}\left(\bm{\upmu},\textbf{p},\textbf{H}\right)+\lambda\phi_{\rm
p}\left(\bm{\upmu}\right)$ (1)

Here, $\Phi$ denotes a general Bayesian cost function with $\phi_{\rm f}$ the
data fidelity term and $\phi_{\rm p}$ the priors. Here, p is the a CT
projection data that is very noisy from a low dose scan, and H the system
matrix of CT, $\bm{\upmu}$ denotes the estimation of linear attenuation map,
i.e. be the output of network. The hyper-paramete $\lambda$ in Eq.1 is
commonly used to adjust the strength of two terms. Hence, with known imaging
system matrix and projection data, the cost function avoids using any ground-
truth label and model imaging physics instead.

In X-ray CT imaging, the data are generally modelled to be bin-to-bin
independent Poisson or Compound Poisson distributed[15]. After negative log
operation of Poisson distributed data, a sinogram projection is approximately
Gaussian with its variance dependent on mean [16]. Thus, we have the
likelihood to be:

$\textbf{P}\left(\textbf{p};\bm{\upmu}\right)=\exp\left(-(\textbf{H}\bm{\upmu}-\textbf{p}\right)^{T}\bm{\Sigma}^{-1}\left(\textbf{H}\bm{\upmu}-\textbf{p})\right)$
(2)

Here, $\bm{\Sigma}^{-1}$ is a diagonal matrix and its diagonal elements denote
the variance of projections. Various priors have been studied in X-ray CT
reconstruction [9, 10, 17, 18, 19]. In this work, we adopt a constraint
derived from none local mean (NLM) [8]. The NLM term can be formulated as:

$\phi_{\rm
NLM}\left(\bm{\upmu}\right)=\frac{1}{2}\Sigma_{i}\Sigma_{i^{\prime}\in\mathcal{N}\left(i\right)}W_{ii^{\prime}}\left(\mu_{i}-\mu_{i^{\prime}}\right)^{2}$
(3)

Here, $\mathcal{N}\left(i\right)$ denotes the neighboring pixels of $i$ within
a searching window and $W_{ii^{\prime}}$ denotes the weight of similarity
computed from:

$W_{ii^{\prime}}=\frac{{\rm
exp}\left(-||\Pi\left(\mu_{i}\right)-\Pi\left(\mu_{i^{\prime}}\right)||_{2,a}^{2}/h\right)}{\Sigma_{j\in\mathcal{N}\left(i\right)}{\rm
exp}\left(-||\Pi\left(\mu_{i}\right)-\Pi\left(\mu_{j}\right)||_{2,a}^{2}/h\right)}$
(4)

With $\Pi\left(\cdot\right)$ being a patch extractor, $||\cdot||_{2,a}$ the
L2-norm with a Gaussian kernel $a$ to weight the distance in a patch and a
normalization parameter. If neighboring patches of two pixels are similar, the
two pixels are likely to be the same kind of materials and gets a high weight.
In this way, NLM can preserve structures well while restrain noise. Combining
Eq. 2 and 3, we obtain a Bayesian cost function which falls into the category
of penalized weighted least square (PWLS). It is used to guide our deep
learning network:

$\Phi\left(\bm{\upmu}\right)=\left(\textbf{p}-\textbf{H}\bm{\upmu}\right)^{T}\bm{\Sigma}^{-1}\left(\textbf{p}-\textbf{H}\bm{\upmu}\right)+\lambda\phi_{\rm
NLM}\left(\bm{\upmu}\right)$ (5)

###  2.2 Network architecture and training

To enforce the cost function as in Eq.5, we constructed a reconstruction
network as displayed in Fig. 1. It takes noisy projection data as inputs and
output reconstructed images at the end. The network is consisted of three
blocks to do dual-domain learning similar to our previous work [20, 21].

Figure 1: Architecture of model-based deep learning network for CT
reconstruction.

First, projection block reduces the noise in projection domain. Projection
block contains five convolutional layers. Each convolutional layer is followed
with leaky-Relu activation function. Next, domain transfer block is an
analytical inverse operator. We adopt the matrix format of FBP reconstruction
as Eq.6

$\bm{\upmu}=\textbf{H}_{W}^{T}\textbf{F}\left(\textbf{W $\cdot$ p}\right)$ (6)

Here, W denotes projection weighting matrix, it weights the projection. F is a
matrix form of RL filtration, it filters the projection data along the
detector axis. $\textbf{H}_{W}^{T}$ denotes the operator of weighted back-
projection. All the matrix in the reconstruction block is pre-calculated and
not trained in this study, but gradients can path through the reconstruction
block from image domain back to projection domain, which is critical in our
network. Following reconstruction block, an image domain U-Net [22] refines
the reconstruction results and output the final estimation. The U-Net
structure achieves big receptive-field with pooling layers, while it recovers
local spatial resolution by combining up-sampled features with un-pooled
features. Both image domain network and projection domain network adopt
residual structure to ease the network training process.

In our training, the network parameters are updated to minimize the Bayesian
cost function over the whole training set.

$v_{1},...,v_{M}=\mathop{\rm argmin}\limits_{v_{1},...,v_{M}}\Sigma_{k\in
1:K}\Phi\left(\bm{\upmu}^{k}\right)$ (7)

Here, $k$ indexes the attenuation map in the training dataset of size $K$, and
$v_{1},...,v_{m},...,v_{M}$ denote the parameters in the network.

For each training step, gradients of cost function with respect to $k^{th}$
output image can be calculated from Eq. 5 as:

$\frac{\partial\Phi\left(\bm{\upmu}^{k}\right)}{\partial\mu_{i}^{k}}=\left[\textbf{H}^{T}\left(\textbf{H}\bm{\upmu}^{k}-\textbf{p}\right)\right]_{i}\bm{\Sigma}_{i}^{-1}+\lambda\Sigma_{i^{\prime}\in\mathcal{N}\left(i\right)}W_{ii^{\prime}}\left(\mu_{i}^{k}-\mu_{i^{\prime}}^{k}\right)$
(8)

Here, subscript $i$ denotes the $i^{th}$ element of a vector. Further, the
gradients of the cost function with respect to network parameters can be
calculated by chain rule:

$\frac{\partial\Phi\left(\bm{\upmu}^{k}\right)}{\partial
v_{m}^{n}}=\Sigma_{i}\frac{\partial\Phi\left(\bm{\upmu}^{k}\right)}{\partial\mu_{i}^{k}}\cdot\frac{\partial\mu_{i}^{k}}{\partial
v_{m}^{n}}$ (9)

with $n$ indexing the $n^{th}$ training iteration.

Parameters are updated with standard gradient descent method to minimize the
cost function:

$v_{m}^{n+1}=v_{m}^{n}-\alpha\Sigma_{k\in
1:K}\frac{\partial\Phi\left(\bm{\upmu}^{k}\right)}{\partial v_{m}^{n}}$ (10)

Here, $\alpha$ denotes the learning rate. Other optimization method can also
be used to replace the gradient descending method. It is worth to emphasize
that we train the network by directly updating parameters of network as above.

Different from supervised pixel-to-pixel loss, Bayesian cost function does not
provide clear optimization direction in image domain for output images. To be
more specific, the gradients of output image do not force output image with a
clear goal in image domain in different training steps, which makes the
network converges slowly in training. Thus, we applied some tricks to
accelerate training process. In the early epochs, we train each data in
training set for several times continuously, then change to another training
data. Training each data continuously help network to catch more robust rule
from Bayesian cost function. Moreover, in the early epochs, we only train the
network with data fidelity term so that the main structure of reconstructed
images can be learned efficiently. After the pre-training process, prior
constrain is added.

##  3 Experiments and results

###  3.1 Study with simulated projection data

We first carried out our experiments on 3872 practical high quality dental CT
images. We simulated the LDCT projections as network inputs. Among them, 3400
randomly chosen simulated projection data were used as training set while the
left projections and the corresponding high-quality images formed test set.
The high-quality dental CT images were not available to training process but
only used as ground-truth for performance evaluation of the reconstruction
results during inference, thus the results can be evaluated both in vision and
quantitatively. Both ground-truth images and reconstruction images were
$640^{2}$ with pixel size 0.25mm. The projections were of uniformly
distributed 720 fan-beam views over 2$\pi$, and 640 detector bins with bin
size 0.512mm. The projections were simulated with ray-tracing forward
projection method and Poisson noise was added with 5000 incident photons.

In Fig.2, ground-truth, LDCT FBP reconstruction, NLM-WLS iterative
reconstruction, and MBDL network reconstruction were displayed for visual
comparison. Details denoted by red boxes were zoomed-in and displayed in the
lower left corner. As shown in Fig.2, the noise was well suppressed in both
iterative reconstruction and MBDL reconstruction, the details of teeth and
bone were restored to a large extent as in zoomed-in.

Figure 2: Comparison of different reconstruction methods. (a): Ground-truth,
(b): LDCT FBP, (c): NLM-WLS iterative reconstruction, (d): LDCT MBDL
reconstruction, (e): Residual map between (a) and (b), (f): Residual map
between (a) and (c), (g): Residual map between (a) and (d).

For further quantitative analysis, we compared different LDCT reconstruction
results with ground-truth images by calculating mean RRMSE and SSIM indexes
based on 25 randomly chosen test images:

${\rm
RRMSE}\left(\bm{\hat{\upmu}},\bm{\upmu}\right)=\frac{||\bm{\hat{\upmu}}-\bm{\upmu}||_{2}}{||\bm{\upmu}||_{2}}$
(11)

${\rm
SSIM}\left(\bm{\hat{\upmu}},\bm{\upmu}\right)=\frac{\left(2\bm{\hat{\upmu}}\bm{\upmu}+C_{1}\right)\left(2\sigma_{\bm{\hat{\upmu}},\bm{\upmu}}^{2}+C_{2}\right)}{\left(\bm{\hat{\upmu}}^{2}+\bm{\upmu}^{2}+C_{1}\right)\left(\sigma_{\bm{\hat{\upmu}}}^{2}+\sigma_{\bm{\upmu}^{2}}^{2}+C_{2}\right)}$
(12)

Here, $\bm{\hat{\upmu}}$ denotes LDCT reconstruction results of different
methods to be evaluated, and $\bm{\upmu}$ corresponding ground-truth.
$C_{1},C_{2}$ are constants. The mean and standard deviation of RRMSE and SSIM
indexes between different FBP reconstruction methods and ground-truth were
compared in Tab.1. Both WLS-NLM iterative reconstruction and MBDL
reconstruction significantly reduce the RRMSE and improve SSIM indexes from
LDCT FBP reconstruction results, while the MBDL reconstructions achieved
slightly better quantitative performance than WLS-NLM iterative
reconstructions in practice.

Table 1:  RRMSE and SSIM between gorund-truth and LDCT reconstruction Method | LDCT FBP | LDCT WLS+NLM IR | LDCT MBDL
---|---|---|---
RRMSE | 0.2560$\pm$0.0094 | 0.1164$\pm$0.0022 | 0.1082$\pm$0.0039
SSIM | 0.9394$\pm$0.0083 | 0.9868$\pm$0.0015 | 0.9885$\pm$ 0.0015

###  3.2 Reconstruction of practical raw data with additional simulated noise

Figure 3: Comparison of different reconstruction methods. (a): NDCT Hanning +
FBP, (b): LDCT Hanning + FBP, (c): NLM-WLS iterative reconstruction, (d): LDCT
MBDL reconstruction.

We further carried out our experiments on practical dental CT raw data while
manually increasing the noise level by additional Poisson noise, thus the
situation is close to reality LDCT. In total, raw data from 3872 patients were
collected. Among them, randomly chosen 3400 data were used as training set
with the others forming test set for the proposed network reconstruction
method. The projections were of uniformly distributed 600 fan-beam views over
2$\pi$, and 656 detector bins covering half of the ROI (reduced size detector)
[23]. The size of detector bins was 0.25mm. To simulate low dose situation,
original normal-dose projection data (no image domain ground truth) were
treated as noise-free reference data, Poisson noise was added with 5000
incident photons. Reconstruction images were $360^{2}$ with pixel size 0.5mm.

In Fig.3, FBP reconstruction of NDCT with Hanning filter, LDCT FBP
reconstruction with Hanning filter, NLM-WLS iterative reconstruction, and MBDL
network reconstruction were displayed for visual comparison. We applied a
5-order Hanning window filtration in FBP reconstructions for noise reduction
in this comparison. Details denoted by red arrow were zoomed-in and displayed
in the lower right corner. Similar with full-simulation cases, both iterative
reconstruction and MBDL network reconstruction suppressed the noise and
restored rich structures of teeth. The two algorithms achieved similar visual
results.

##  4 Conclusion

In this work, we proposed a model-based deep learning reconstruction for LDCT.
The proposed method adopt a Bayesian cost function and construct a dual-domain
network allowing easy error propagation for the cost function minimization.
The MBDL network is trained to minimize the cost function for the whole
dataset with no ground-truth information needed. After training, the network
can complete reconstruction without time-consuming iterative steps and reduced
noise in a level comparable to iterative methods. Fundamentally, the proposed
method integrates the optimization principle of iterative reconstruction
methods with deep learning. This method can be easily extended to other scan
modality and geometry with constraints or various combination of
regularizations for problems such as sparse-view and limited-angle CT
reconstruction.

## References

  * [1] Ej Hall and Dj Brenner.  Cancer risks from diagnostic radiology.  British Journal of Radiology, 81(965):362–378, 2008.
  * [2] D P Naidich, C H Marshall, . Gribbin, C., R S Arams, and D I Mccauley.  Low-dose ct of the lungs: preliminary observations.  Radiology, 175(3):729, 1990.
  * [3] Funama Yoshinori, Awai Kazuo, Miyazaki Osamu, Nakayama Yoshiharu, Goto Taiga, Omi Yasuo, Shimonobo Toshiaki, Liu Duo, Yamashita Yasuyuki, and Hori Shinichi.  Improvement of low-contrast detectability in low-dose hepatic multidetector computed tomography using a novel adaptive filter: evaluation with a computer-simulated liver including tumors.  Investigative Radiology, 41(1):1–7, 2006.
  * [4] Funama Yoshinori, Awai Kazuo, Nakayama Yoshiharu, and et.al.  Radiation dose reduction without degradation of low-contrast detectability at abdominal multisection ct with a low–tube voltage technique: phantom study.  Radiology, 237:905–910, 2005.
  * [5] Wang Jing, Hongbing Lu, Tianfang Li, and Zhengrong Liang.  Sinogram noise reduction for low-dose ct by statistics-based nonlinear filters.  In Medical Imaging: Image Processing, 2005.
  * [6] Anja Borsdorf, Steffen Kappler, Rainer Raupach, Frederic Noo, and Joachim Hornegger.  Local orientation-dependent noise propagation for anisotropic denoising of ct-images.  In Nuclear Science Symposium Conference Record, 2010.
  * [7] Y. Zhang, J. Rong, H. Lu, Y. Xing, and J. Meng.  Low-dose lung ct image restoration using adaptive prior features from full-dose training database.  IEEE Trans Med Imaging, 36(12):2510–2623, 2017.
  * [8] Zhang Hao, Jianhua Ma, Wang Jing, Liu Yan, and Zhengrong Liang.  Adaptive nonlocal means-based regularization for statistical image reconstruction of low-dose x-ray ct.  In Medical Imaging: Physics of Medical Imaging, 2015.
  * [9] Hao Zhang, Jianhua Ma, Jing Wang, Yan Liu, Hongbing Lu, and Zhengrong Liang.  Statistical image reconstruction for low-dose ct using nonlocal means-based regularization.  Computerized Medical Imaging & Graphics, 38(6):423–435, 2014.
  * [10] Xu Qiong, Yu Hengyong, Mou Xuanqin, Zhang Lei, Hsieh Jiang, and Wang Ge.  Low-dose x-ray ct reconstruction via dictionary learning.  IEEE Transactions on Medical Imaging, 31(9):1682–1697, 2012.
  * [11] Ge Wang, Hu Chen, Jiliu Zhou, Ke Li, Peixi Liao, Weihua Zhang, and Yi Zhang.  Low-dose ct via convolutional neural network.  Biomedical Optics Express, 8(2):679, 2017.
  * [12] H. Chen, Y. Zhang, M. K. Kalra, F. Lin, Y. Chen, P. Liao, J. Zhou, and G. Wang.  Low-dose ct with a residual encoder-decoder convolutional neural network (red-cnn).  IEEE Trans Med Imaging, PP(99):2524–2535, 2017.
  * [13] C. Zhang and Y. Xing.  Ct artifact reduction via u-net cnn.  In Society of Photo-optical Instrumentation Engineers, 2018.
  * [14] Q. Yang, P. Yan, Y. Zhang, H. Yu, Y. Shi, X. Mou, M. K. Kalra, Y. Zhang, L. Sun, and G. Wang.  Low-dose ct image denoising using a generative adversarial network with wasserstein distance and perceptual loss.  IEEE Transactions on Medical Imaging, 37(6):1348–1357, 2017.
  * [15] Idris A. Elbakri and Jeffrey A. Fessler.  Efficient and accurate likelihood for iterative image reconstruction in x-ray computed tomography.  Proc Spie, 5032, 2003.
  * [16] J Wang, Z Lu, Hliang, D Eremina, G Zhang, S Wang, J Chen, and J Manzione.  An experimental study on the noise properties of x-ray ct sinogram data in radon space.  Physics in Medicine & Biology, 53(12):3327–3341, 2008.
  * [17] Ritschl Ludwig, Bergner Frank, Fleischmann Christof, and Kachelriess Marc.  Improved total variation-based ct image reconstruction applied to clinical data.  Physics in Medicine & Biology, 56(6):1545–1561, 2011.
  * [18] Kim Kyungsang, Ye Jong Chul, Worstell William, Ouyang Jinsong, Rakvongthai Yothin, El Fakhri Georges, and Li Quanzheng.  Sparse-view spectral ct reconstruction using spectral patch-based low-rank penalty.  IEEE Transactions on Medical Imaging, 34(3):748–760, 2015.
  * [19] Ruoqiao Zhang, Hye Ye Dong, Debashish Pal, Jean Baptiste Thibault, and Charles A. Bouman.  A gaussian mixture mrf for model-based iterative reconstruction with applications to low-dose x-ray ct.  IEEE Transactions on Computational Imaging, 2(3):359–374, 2016\.
  * [20] Kaichao Liang, Hongkai Yang, and Yuxiang Xing.  Comparison of projection domain, image domain, and comprehensive deep learning for sparse-view x-ray ct image reconstruction.  Arxiv 1804.04289, 2018.
  * [21] Hong Kai Yang, Kai Chao Liang, Ke Jun Kang, and Yu Xiang Xing.  Slice-wise reconstruction for low-dose cone-beam ct using a deep residual convolutional neural network.  Nuclear Science and Engineering, 30(4):53–61, 2019.
  * [22] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In International Conference on Medical Image Computing & Computer-assisted Intervention, 2015.
  * [23] Liang Li, Zhiqiang Chen, Li Zhang, Yuxiang Xing, and Kejun Kang.  A cone-beam tomography system with a reduced size planar detector: A backprojection-filtration reconstruction algorithm as well as numerical and practical experiments.  Applied Radiation & Isotopes Including Data Instrumentation & Methods for Use in Agriculture Industry & Medicine, 65(9):1041–1047, 2007.