# Benchmark Results for Bookshelf Organization Problem as Mixed Integer
Nonlinear Program with Mode Switch and Collision Avoidance

Xuan Lin1, Gabriel I. Fernandez1, and Dennis W. Hong1 1X. Lin, Gabriel I.
Fernandez, and Dennis W. Hong are with the Robotics and Mechanisms Laboratory,
Department of Mechanical and Aerospace Engineering, University of California,
Los Angeles, CA 90095, USA<EMAIL_ADDRESS>
###### Abstract

Mixed integer convex and nonlinear programs, MICP and MINLP, are expressive
but require long solving times. Recent work that combines data-driven methods
on solver heuristics has shown potential to overcome this issue allowing for
applications on larger scale practical problems. To solve mixed-integer
bilinear programs online with data-driven methods, several formulations exist
including mathematical programming with complementary constraints (MPCC),
mixed-integer programming (MIP). In this work, we benchmark the performances
of those data-driven schemes on a bookshelf organization problem that has
discrete mode switch and collision avoidance constraints. The success rate,
optimal cost and solving time are compared along with non-data-driven methods.
Our proposed methods are demonstrated as a high level planner for a robotic
arm for the bookshelf problem.

##  I Introduction

Optimization-based methods are useful tools for solving robotic motion
planning problems. Typical approaches such as mixed-integer convex programs
(MICPs) [1, 2], nonlinear or nonconvex programs (NLPs) [3, 4, 5] and mixed-
integer NLPs (MINLPs) [6] offer powerful tools to formulate these problems.
However, each has its own drawbacks. NLPs tend to converge to local optimal
solutions. In practice, local optimal solutions can have bad properties, such
as inconsistent behavior as they depend on initial guesses. Mixed-integer
programs (MIPs) explicitly treat the discrete variables. Branch-and-bound is
usually used to solve MIPs [7]. MIP solvers seek global optimal solutions with
more consistent behavior than NLP solvers. For small-scale problems, these
algorithms usually find optimal solutions within a reasonable time [2, 8]. On
the contrary, MIPs can require impractically long solving times for problems
with a large number of integer variables [9]. MINLPs incorporate both integer
variables and nonlinear constraints, hence, very expressive. Unfortunately, we
lack efficient algorithms to tackle MINLPs. As a result, it is difficult to
implement most of optimization schemes online.

Figure 1: Angle and x position of solutions to the bookshelf problem’s
furthest left book are depicted by colored numbers indicating cluster
membership. Each box on the grid corresponds to a different relaxation scheme
using integer variables: (a) full relaxation of the whole space, (b)
relaxation over the space with data, (c) relaxation for specific modes from
clusters. For visualization purposes only mode 4 region is explicitly
displayed. Notice how different modes can reduce the amount of integer
variables needed.

Recently, researchers have started to investigate machine learning methods to
gather problem specific heuristics and speed up the optimization solving
process. Standard algorithms to solve MIPs such as branch-and-bound and
cutting plane methods rely on heuristics to quickly remove infeasible regions.
Learning methods can be used to acquire better heuristics. For example, [10]
used graph neural networks to learn heuristics. [11] used reinforcement
learning to discover efficient cutting planes. On the other hand, data can be
collected to learn and solve specific problems [12, 13]. In [13], the authors
proposed CoCo which collects problem inputs and solved integer variable
strategies offline and trains a neural network. Effectively, this becomes a
classification problem, where each strategy has a unique label. For online
solving, the neural network will propose candidate solutions, reducing it to a
convex program. The results show that CoCo can solve MIPs with around 50
integer variables within a second. However, CoCo has relatively simple neural-
network structure. Since the number of possible strategies equals to
$2^{\lambda}$, where $\lambda$ is the number of integer variables, the number
of unique strategies tend to be close to the total amount of data. This
induces overfitting.

Collision avoidance constraints typically appear in multi-agent motion
planning problems. Single-agent collision avoidance with the stationary
environment is relatively straightforward as the collision-free regions can be
precomputed such that the moving agent stays within them. Polyhedrons are
typical models for collision-free environment [14]. If moving agents have
simple shapes such as circles, collision-free constraints can simply enforce
that distance between two agents are larger than a known value [15]. However,
if moving agents have complex shapes, collision avoidance becomes more
challenging. Naive methods include enforcing collision avoidance between a
large number of points sampled within each agent. This approach is obviously
computationally expensive. Methods based on KKT conditions are proposed [16]
to enforce collision avoidance between convex shapes. This method does not
allow the agent to make contact without sacrificing precision. In this paper,
we propose a novel approach using separating planes to enforce collision
avoidance between convex agents but still allow them to make precise contacts.
If the agent has a non-convex shape, it can be decomposed into convex shapes
such that this approach can still be used.

For many multi-agent optimization problems, agents can operate under various
modes encoded with discrete variables. The discrete mode switching with
collision avoidance constraints leads to a mixed-integer bilinear formulation
which is challenging to solve. Leveraging on pre-solved data is one approach
to solve these formulations fast and reliably. To get a comprehensive
understanding of the performance of various data-driven methods on such
problems, we introduce the bookshelf organization problem for a benchmark.
Given a bookshelf with several books on top, an additional book needs to be
placed on the shelf with minimal disturbance to the existing books. The
bookshelf problem works well as a good benchmark because it: 1) is an MINLP
that can be converted to an MICP problem with hundreds of integer variables,
2) can easily be scaled to push algorithms to their limit, and 3) has
practical significance where data can reasonably be gathered, such as in the
logistics industry.

To summarize, our contributions are as follows:

  1. 1.

Formulate the bookshelf organization problem as an MINLP and solve it within
seconds.

  2. 2.

Benchmark the performance of different data-driven formulations on the
bookshelf organization problem.

  3. 3.

Implement the solutions on hardware to stow items into the shelf.

##  II Related Works

###  II-A Motion Library

Motion library is a typical method to leverage offline computing power. A
large number of controllers are pre-computed offline and stored inside a
library. Online, the controllers are selected based on the robot state in the
nearest neighbor manner. This method can provide more robust and global
policies faster than online solving techniques such as dynamic programming
[17]. Motion libraries have been used to synthesize controllers for bipedal
walking [18], and simulated multi-link swimming [19]. However, this method
does not adapt well to situations not previously seen in the library.

###  II-B Parametric Programming

Parametric programming is a technique to build a function that maps to the
optimization solutions from varying parameters of the optimization problem
[20, 21]. Previous works has investigated parametric programming on linear
programming [22], quadratic programming [23], mixed-integer nonlinear
programming [24]. As an implementation of parametric programming on controller
design, explicit MPC [23, 25] tries to solve the model-predictive control
problem offline, and stores active sets. When computed online, the problem
parameter can be used to retrieve the active sets. [23] solved a constrained
linear quadratic regulator problem with explicit MPC and proved that the
active sets are polyhedrons. Therefore, the parameter space can be partitioned
using an algorithm proposed by [26]. However, when the problem is non-convex,
the active sets do not in general form polyhedrons. Therefore, it is
significantly more complex to compute critical regions offline.

###  II-C Learning Problem-solution Mapping

The basic ideas behind parametric optimization is to construct a mapping from
the descriptor of the problem to the solutions of the problem. For simple
problems, this may be done analytically with parametric optimization
techiques. For more challenging problems, it makes more sense to use learning
approaches. For relatively smaller datasets, previous works [27, 28, 12]
directly store the data points and pick out warm-start using non-parametric
learning such as K-nearest neighbor (KNN), and solve the online formulation.
Effectively, this approach adds an online adaptation step to the motion
library method, which turns out to work well. On the other hand, modern
learning techniques such as neural-networks can learn an embedding of a larger
set of parameters that maps to the solutions [13]. One advantage of using
neural-network methods is the capability to deal with out-of-distribution
situations that are not included in the training set [29].

###  II-D Mixed-integer Programming

To allow online adaptation, one needs to solve an online optimization
formulation in real-time. This online formulation can have several forms.
Typical robotics problem includes discrete contact and discrete operating
modes. There are generally two ways to include discrete variables in the
optimization formulation. Mixed-integer programs include continuous variables
and explicit discrete variables. The definition of mixed-integer variables is
that if the discrete variables are relaxed into continuous ones, the problem
becomes convex [30]. Typical mixed-integer programming solvers use branch-and-
bound methods [7], cutting plane methods [31], and rely on heuristics to
choose the branch to explore. Despite the worst case of solving time, a large
portion of problems only requires exploring a small portion of the search tree
[32]. Mixed-integer programs have been implemented for online motion planning
such as [8].

###  II-E Mathematical Programming with Complementary Constraints

On the other hand, mathematical programs with complementary constraints (MPCC)
models discrete modes through continuous variables with complementary
constraints. Complementary constraints enforce a pair of variables such that
if one of them is non-zero, the other one should be zero. This constraint is
traditionally hard to solve and is sensitive to initial guesses. Algorithms
such as time-stepping [33], pivoting [34], central path methods [35] are
proposed to resolve complementarity. In the robotics community, complementary
constraints are typically used to optimize over gaits for trajectory
optimization [36, 37] or control with implicit contacts [38] where contact
forces and distance to the ground are complementary with each other. On the
other hand, complementary constraints can also be used to model binary
variables [39].

##  III Bookshelf Organization Problem Setup

Figure 2: Complete formulation of the bookshelf organization problem.

Assume a 2D bookshelf with limited width $W$ and height $H$ contains
rectangular books where book $i$ has width $W_{i}$ and height $H_{i}$ for
$i=1,...,N-1$. A new book, $i=N$, is to be inserted into the shelf. The
bookshelf contains enough books in various orientations, where in order to
insert book $N$, the other $N-1$ books may need to be moved, i.e., optimize
for minimal movement of $N-1$ books. This paper focuses on inserting one book
using a single robot motion during which the book is always grasped. The
problem settings can be extended to a sequence of motions with re-grasping.
This problem is useful in the logistics industry.

Fig. 2 shows the constraints, variables and objective function for the
bookshelf problem. The variables that characterize book $i$ are: position
$\textbf{x}_{i}=[x_{i},y_{i}]$ and angle $\theta_{i}$ about its centroid.
$\theta_{i}=0$ when a book stands upright. The rotation matrix is:
$\textbf{R}_{i}=[cos(\theta_{i}),\ -sin(\theta_{i});\ sin(\theta_{i}),\
cos(\theta_{i})]$. Let the 4 vertices of book $i$ be $\textbf{v}_{i,k}$,
$k=1,2,3,4$. The constraint A in Fig. 2 shows the linear relationship between
$\textbf{x}_{i}$ and $\textbf{v}_{i,k}$, where $\textbf{h}_{i,k}$ is the
constant offset vector from its centroid to vertices. Constraint B enforces
that all vertices of all books stay within the bookshelf, a linear constraint.
Constraint C enforces the orthogonality of the rotation matrix, a bilinear
(non-convex) constraint. Constraint D enforces that the angle $\theta_{i}$
stays within $[-90\degree,90\degree]$, storing books right side up.

To ensure that the final book positions and orientations do not overlap with
each other, separating plane constraints are enforced. For convex shapes, the
two shapes do not overlap with each other if and only if there exists a
separating hyperplane $\textbf{a}^{T}\textbf{x}=b$ in between [40]. That is,
for any point $\textbf{p}_{1}$ inside shape 1 then
$\textbf{a}^{T}\textbf{p}_{1}\leq b$, and for any point $\textbf{p}_{2}$
inside shape 2 then $\textbf{a}^{T}\textbf{p}_{2}\geq b$. This is represented
by constraint E. Constraint F enforces a to be a normal vector. Both E and F
are bilinear constraints.

Finally, we need to assign a state to each book $i$. For each book, it can be
standing straight up, laying down on its left or right, or leaning towards
left or right against some other book, as shown in the far right column in
Fig. 2. For each book $i$, we assign a set of integer variables ${z_{i}}$. If
book $i$ stands upright ($z_{i,2}=1$) or lays flat on its left ($z_{i,1}=1$)
or right ($z_{i,3}=1$), constraints I1 or J1 or H1 are enforced, respectively.
If book $i$ leans against another book (or left/right wall which can be
treated as static books) on the left or right, constraints in K and L are
enforced, respectively. To this end checks need to indicate the contact
between books. By looking at the right column in Fig. 2, we can reasonably
assume that the separating plane $\textbf{a}^{T}\textbf{x}=b$ always crosses
vertex 1 of the book on the left and vertex 4 of the book on the right. This
is represented by bilinear constraints, K1 and L1. In addition, the books need
to remain stable given gravity. Constraints K2 and L2 enforce that a book is
stable if its $x$ position stays between the supporting point of itself
(vertex 2 if leaning rightward and vertex 3 if leaning leftward) and the $x$
position of the book that it is leaning onto. Lastly, constraint K3 and L3
enforce that the books have contact with the _ground_. For practical reasons,
we assume that books cannot stack onto each other, i.e, each book has to touch
the _ground_ of bookshelf at at least one point. We note that constraints in
H, I, J, K, and L can be easily formulated as MICP constraints using big-M
formulation [41], such that they are enforced only if the associated integer
variable $z=1$. Also, it can easily be extended to allow stacking for our
problem. Any contact conditions between pairs of books may also be added into
this problem as long as it can be formulated as mixed-integer convex
constraint. Overall, this is a problem with integer variables, ${z_{i}}$, and
non-convex constraints C, E, F, K1, and L1, hence, an MINLP problem.

Practically, this problem presents challenges for retrieving high quality
solutions. If robots were used to store books, the permissible solving time is
several seconds, and less optimal solutions means longer realization times.
The video of our hardware implementation demonstrates non-optimal insertion
cases that requires multiple additional robot motions that dramatically
increase the chance of failure. There are several potential approaches to
resolving this issue: fix one set of nonlinear variables and solve MICP [42],
convert the nonlinear constraints into piece-wise linear constraints and
formulate them into an MICP [1], or directly apply MINLP solvers such as
BONMIN [43]. As expected, these approaches struggle to meet the requirements.

##  IV Learning Algorithm

Assume that we are given a set of problems parametrized by $\Theta$ that is
drawn from a distribution $D(\Theta)$. For each $\Theta$, we seek a solution
to the optimization problem:

$\displaystyle\underset{\textbf{x},\ \textbf{z}}{\text{minimize}}\
f_{obj}(\textbf{x},\textbf{z};\Theta)$ (1) s. t. $\displaystyle
f_{i}(\textbf{x},\textbf{z};\Theta)\leq 0,\ \ i=1,...,m_{f}$ $\displaystyle
b_{j}(\textbf{x},\textbf{z};\Theta)\leq 0,\ \ j=1,...,m_{b}$

Where x denotes continuous variables and z binary variables with
$z_{i}\in\\{0,1\\}$ for $i=1,...,dim(\textbf{z})$. Constraints $f_{i}$ are
mixed-integer convex, meaning if the binary variables z are relaxed into
continuous variables $\textbf{z}\in[0,1]$, $f_{i}$ becomes convex. Constraints
$b_{j}$ are mixed-integer bilinear, meaning that relaxing the binary variables
gives bilinear constraints. Without loss of generality, x and z are assumed to
be involved in each constraint. We omit equality constraints in (1) as they
can be turned into two inequality constraints from opposite directions.

In general, there are two directions to convert a mixed-integer bilinear
program to either a mixed-integer linear program or a nonlinear program. To
convert to a mixed-integer linear program, one needs to convert non-convex
constraints to mixed-integer linear constraints. For example, we can convert
trigonometry constraints into piece-wise linear constraints [1], or convert
bilinear constraints into mixed-integer envelope constraints [30]. However,
mixed-integer envelope constraint tend to significantly increase the solving
time. In [9], the trajectory planner is based on mixed-integer formulation
with envelopes to approximate the bilinear moment constraints. To walk a few
steps upstairs, the solver can take several hours to find a feasible solution.
The other approach is to convert the binary variables into continuous
variables with complementary constraints. This approach requires the solver to
deal with a potentially large number of complementary constraints known to be
computationally difficult. Directly solving the NLP with complementary
constraints without an informed initial guess result in high chance of
infeasibility, unless special treatment such as [33] are used. Those methods,
however, result in slower solving speed from our benchmark results.

For both directions mentioned above, incorporating pre-solved data to learn a
good warm-start online can improve the solving speed or rate of feasibility,
as we will show in our experiments. If the learning agent with the problem
description $\Theta$ can partially provide feasible integer variables, the
size of the mixed-integer optimization problem can be reduced and hence solved
faster. If the learning agent can provide a complete list of integer
variables, the problem reduces to a convex optimization permitting convex
solvers such as OSQP [44]. On the other hand, if the learning agent provides a
good initial guess for the variables involved in complementary constraints,
the feasibility to solve the problem significantly increases.

To better understand the performances of the approaches mentioned above, we
compared their performances on the bookshelf problem. In addition, we have
also tested the performance of non-data-driven methods, including manually
designed heuristics, and nonlinear ADMM [45] as benchmarks. We first collect
the initial dataset. This is done by sampling the problem parameters $\Theta$,
then solving the samples with a non-data-driven approach. Either manually
designed heuristics or nonlinear ADMM can be used. We can then use the initial
dataset as warm-start points to solve the new problem instances and enlarge
the dataset. With warm-start points, data-driven methods can be used, and the
new dataset has better optimal costs. This procedure iterates until we arrive
at the final dataset. The final dataset is used to learn good problem-solution
mapping for warm-start online. The online formulation can be MPCC, MIP, or a
convex formulation. In the following sections, we describe the implementation
details for each online formulation.

TABLE I: Benchmark Results

| ADMM - Knitro | ADMM - IPOPT | MIP - 50 | MIP - 20 |  | MPCC_default
---
\- Knitro
| MPCC_default
---
\- IPOPT
| MPCC_Manual
---
\- Knitro
| MPCC_Manual
---
\- IPOPT
| MPCC_KNN3
---
\- Knitro (15000data)
| MPCC_KNN3
---
\- IPOPT (15000 data)
| MPCC_KNN3 (Rev)
---
\- Knitro (36000 data)
| Convex - 200
---
(15000 data)
| Convex - 600
---
(15000 data)
| Convex - 200
---
(78000 data)
Note |  |  |  | Each cluster
---
50 samples
| Each cluster
---
20 samples
| Using default initial guess from
---
CasADi
Using manually selected 1 initial guess | KNN selects top 3 guesses from 15000 data |  | Same to left but select
---
the worst 3 guesses
| Max iteration
---
of OSQP=200
| Max iteration
---
of OSQP=600
| Currently KNN
---
learner
Success Rate | 96.5% | 94.75% | 100 % | 98.75% | 1.25 % | 0 % | 78.25% | 78.25% | 99.25% | 98.5% | 92.5% | 94.5% | 98.75% | 96%
Avg. Solve Time | 260 ms | 522 ms | 790 ms | 616 ms | 72 ms |  | 29 ms | 81 ms | 30 ms | 96 ms | 61 ms | 65 ms | 80 ms | 16 ms
Max. Solve Time | 1.29 sec | 3.59 sec | 17.3 s | 70 s | 82 ms |  | 198 ms | 1.96 sec | 280 ms | 990 ms | 327 ms | 190 ms | 400 ms | 174 ms
Avg. Objective | -7250 | -7339 | -8606 | -8435 | -5703 |  | -8620 | -8952 | -8258 | -8618 | -6461 | -8631 | -8687 | -8637
Avg. # of iterations | 4.73 | 4.6 | 1 | 1 |  |  | 1 | 1 | 1.05 | 1.06 | 1.25 | 4.85 | 3.85 | 2.57
Solver | Gurobi+Knitro | Gurobi+IPOPT | Gurobi | Gurobi | Knitro | IPOPT | Knitro | IPOPT | Knitro | IPOPT | Knitro | OSQP | OSQP | OSQP

###  IV-A Complementary Formulation

An equivalent formulation of the bookshelf organization problem formulation is
to turn all the integer variables $z_{i}$ into continuous variables with
complementary constraints, such that the complete formulation can be solved
through NLP solvers. We implemented the formulation from [39]:

$z_{i}(1-z_{i})=0$ (2)

Which is equivalent to the constraint $z_{i}\in\\{0,1\\}$. We have also tried
other implementations in [46] and get similar or worse performances.

After collecting the dataset, we can learn the optimal problem-solution
mapping from the problem parameter $\Theta$ to the solution
$(\textbf{x},\textbf{z})$. Online, the learner samples multiple warm-start
points $(\textbf{x},\textbf{z})$ for the NLP. The NLP with complementary
constraints are then solved using the sampled warm-start one-by-one until a
feasible solution is returned.

###  IV-B MIP Formulation

We convert bilinear constraints into mixed-integer linear constraints by
gridding the solution space $\mathbb{S}$ and approximating the constraints
locally inside grids with McCormick envelopes similar to [30]. A McCormick
envelope relaxation of one bilinear constraint $w=xy$ [47] is the best linear
approximation defined over a pair of lower and upper bounds $[x^{L},x^{U}]$
and $[y^{L},y^{U}]$. Therefore, we first assign grids $G(\textbf{x})$ to
$\mathbb{S}$. Let $\\{g_{l}\\}$, $l=1,...,L$ be one cell in the grid with
upper and lower bounds $[\textbf{x}^{L},\textbf{x}^{U}]$. We introduce
additional integer variables n, $n_{i}\in\\{0,1\\}$. Each unique value of n
corresponds to one cell in the grid within which McCormick envelope
relaxations are applied. The constraints
$b_{j}(\textbf{x},\textbf{z};\Theta)\leq 0$, $j=1,...,m_{b}$ are converted
into $L$ constraints:
$E_{b_{j},l}(\textbf{x},\textbf{z},\textbf{n};\Theta)\leq 0$ for
$j=1,...,m_{b}$ and $l=1,...,L$, turning (1) into an MICP:

$\begin{split}&\underset{\textbf{x},\ \textbf{z},\
\textbf{n}}{\text{minimize}}\
f_{obj}(\textbf{x},\textbf{z},\textbf{n};\Theta)\\\ \text{s.
t.}\qquad\;f_{i}(\textbf{x},\textbf{z};\Theta)\leq 0,&\ \ i=1,...,m_{f}\\\
E_{b_{j},l}(\textbf{x},\textbf{z},\textbf{n};\Theta)&\leq 0,\ \
j=1,...,m_{b},\ \ l=1,...,L\qquad\end{split}$ (3)

We use a $log_{2}N$ formulation which means $N$ grids are represented by
$\lceil log_{2}N\rceil$ binary variables. For example, $17\sim 32$ grids are
represented by 5 integer variables. Since the intervals generated by the
clustering methods are not necessarily connected, the proposed method in [30]
does not work. The formulation we use can incorporate intervals that are not
connected hence more general. The price to pay is additional continuous
variables. Please see Appendix for details of the formulation. Despite the
relatively smaller number of binary variables due to $log_{2}N$ formulation,
approximating nonlinear constraints with mixed-integer convex constraints
still generates a large number of integer variables when high approximation
accuracy is desired. As a result, the formulation suffers from extended
solving time. We use data-driven methods to reduce the number of integer
variables and learn good warm-starts for online applications.

####  IV-B1 Reduced Scale MIP Formulation

[48] proposes ReDUCE as a method to generate smaller scale MIP problems. The
basic idea of using unsupervised learning to reduce the size of mixed-integer
programs is to identify and cluster the important regions on the solution
manifold, and generate MIP formulations with smaller scale. If
$\mathbb{R}^{dim(\textbf{x})}$ is segmented into smaller regions, e.g.
clusters, the required integer variables for each cluster can be reduced. This
may be seen in Fig. 1 which is an instance of 2-dimensional (dim) $X(\Theta)$
from a bookshelf experiment described in Sec. V-A.

This paragraph reviews the general steps of ReDUCE. To begin, ReDUCE uses pre-
solved dataset $(\Theta_{k},\textbf{x}_{k})$, $k=1,...,K$. We pre-assign grids
$G(\textbf{x})$ to the solution space $\mathbb{S}$. The size of grids depends
on the approximation accuracy requirement for bilinear constraints. ReDUCE
begins by performing unsupervised learning on the pre-solved dataset to
retrieve clusters that indicates regions on $X(\theta)$. Density-Based Spatial
Clustering of Applications with Noise (DBSCAN) [49] was used to cluster on x
which gives clusters ${\textbf{x}_{1}},{\textbf{x}_{2}},...,{\textbf{x}_{C}}$,
where $C$ is the number of clusters. We then trace
$\Theta\rightarrow\textbf{x}$ map backwards to create clusters in $\Theta$
space, i.e.
$(\\{\Theta\\}_{1},\\{\textbf{x}\\}_{1}),...,(\\{\Theta\\}_{C},\\{\textbf{x}\\}_{C})$.
Next, a supervised classifier is trained to classify $(\Theta,c)$. We use
Random Forest which requires relatively smaller amounts of data to train than
deep learning methods. For all $\\{\textbf{x}\\}_{c}$ in cluster $c$, we find
the grid cells that they occupy and re-assign integer variables
$\textbf{z}_{c}$ and $\textbf{n}_{c}$ to each cluster. At this point, the
original MIP problem is segmented into smaller scale MIP problems with faster
solving speed. When a new problem instance comes online, the classifier is
used to point the new problem to one of the existing clusters and the reduced
scale MIP is solved.

The idea of solving smaller sub-MIPs bears similarities with algorithms such
as RENS [50] and Neural Diving [10] where a subset of integer variables are
fixed from linear program solutions or a learned model, while the others are
solved. Our method, however, segments the problem through an unsupervised
clustering approach.

####  IV-B2 Convex Formulation

In the extreme case, if the learner can provide a complete warm-start of the
binary variables, the problem online becomes convex permitting fast QP
solvers. [51, 13, 52] define integer strategies to be tuples of
$\mathcal{I}(\Theta_{i})=(\textbf{z}^{*},\textbf{n}^{*},\mathcal{T}(\Theta_{i}))$,
where $(\textbf{x}^{*},\textbf{z}^{*},\textbf{n}^{*})$ is an optimizer for
problem (3), $\mathcal{T}(\Theta_{i})=\\{i\in{1,...,m_{f}},\ l\in
1,...,L|f_{i}(\textbf{x}^{*},\textbf{z}^{*};\Theta)=0,\
E_{b_{j},l}(\textbf{x}^{*},\textbf{z}^{*},\textbf{n}^{*};\Theta)=0\\}$ is the
set of active inequality constraints. Given an optimal integer strategy
$\mathcal{I}(\Theta_{i})$, solutions to (3) can be retrieved through solving a
convex optimization problem. This idea is adopted for learning warm-start with
a convex formulation online. We have implemented both K-nearest neighbor to
sample the points that are closest to the new problem $\Theta$. We have also
tested a simple neural-network structure identical to [13].

###  IV-C ADMM Formulation

We implemented a non-data-driven benchmark based on the ADMM algorithm in
[53]. The idea is to separate the MINLP into a MIP formulation and an NLP
formulation with the exact copies of variables but different constraints, and
iterate between them until the consensus is reached. The MIP formulation
contains constraints $A,B,D,G,H1,I1,J1,K2,K3,L2,L3$, while missing bilinear
constraints $C,E,F,K1,L1$. The NLP formulation contains constraints
$A,B,C,D,E,F,H,I,J,K,L$, while missing the mode constraint G which says
$z_{i}$s are integers. The formulation is elaborated in Eqn. 4.

$\displaystyle\underset{\textbf{x}_{i},\ \textbf{R}_{i},\ \textbf{v}_{i,k},\
\textbf{a}_{j},\ b_{j},\ z_{i}}{\text{minimize}}\quad f_{obj}$ (4) subject to
Mixed integer constraints: A, B, D, G, H1, I1, J1, K2, K3, L2, L3 Nonlinear
constraints: A, B, C, D, E, F, H, I, J, K, L

The ADMM iterates the following 3 steps: first it solves the MIP formulation
in (4), then solves the NLP formulation in (4). Finally it updates the dual
variables. The process continues till convergence. Let $\textbf{var}_{1}$ and
$\textbf{var}_{2}$ be the copies of the variables in formulation (4). Let w be
the dual variables. G are the weights. The first step is:

$\displaystyle\underset{\textbf{var}_{1}}{\text{minimize}}\quad||\textbf{var}_{1}^{i}-\textbf{var}_{2}^{i}+\textbf{w}^{i}||_{\textbf{G}_{k}}$
s.t. Mixed-integer constraints in (4)

The next step solves the NLP:

$\displaystyle\underset{\textbf{var}_{2}}{\text{minimize}}\quad||\textbf{var}_{2}^{i}-(\textbf{var}_{1}^{i+1}+\textbf{w}^{i}))||_{\textbf{G}_{k}}$
s.t. Nonlinear constraints in (4)

Finally, we update the dual variables and weights:

$\displaystyle\textbf{w}^{i+1}=\textbf{w}^{i}+\textbf{var}_{1}^{i}-\textbf{var}_{2}^{i}$
(5) $\displaystyle\textbf{G}_{k+1}=\gamma\textbf{G}_{k}$
$\displaystyle\textbf{w}^{i+1}=\textbf{w}^{i}/\gamma$

##  V Experiment

###  V-A Setup

We place 3 books inside the shelf where 1 additional book is to be inserted.
Grids are assigned to the variables involved in the non-convex constraint C,
E, F, K1, L1: $\textbf{R}_{i}(\theta_{i})$, $\textbf{a}_{j}$ and
$\textbf{v}_{i,k}$. These variables span a 48-dim space. The rotation angles
$\theta_{i}$, which includes $\textbf{R}_{i}$, are gridded at $\frac{\pi}{8}$
intervals. Elements in $\textbf{a}_{j}$ are gridded on 0.25 intervals.
Elements in $\textbf{v}_{i,k}$ are gridded at intervals $\frac{1}{4}$ the
shelf width $W$ and height $H$. Table II lists the number of grids for each
non-convex variable. Our MICP formulation results in 130 integer variables in
total. The input vector includes the center positions, angles, heights and
widths of stored books and height and width of the book to be inserted. The
input dimension is 17.

The problem instances are generated using a 2-dim simulated environment of
books on a shelf. Initially, 4 randomly sized books are arbitrarily placed on
the shelf, and then 1 is randomly removed and regarded as the book to be
inserted. Contrary to the sequence, the initial state with 4 books represents
one feasible (not necessarily optimal) solution to the problem of placing a
book on a shelf with 3 existing books. This guarantees that all problem
instances are feasible. Since this problem can be viewed as high-level
planning for robotic systems, the simulated data is sufficient. For
applications outside of the scope of this paper real-world data may be
preferable in this pipeline.

All methods are tested on 400 randomly sampled bookshelf problems in the same
distribution using the same method mentioned above. All results are listed in
table I.

###  V-B Complementary Formulation

As mentioned in the Section V-B, we turn the binary variables into continuous
variables with complementary constraints. We selected both IPOPT and Knitro
solvers to solve the MPCC formulation. A small relaxation of
$\epsilon=10^{-3}$ on the right hand side of the complementary constraints is
used to increase the chance of getting a feasible solution.

Zero Warm Start As a baseline, we solve the problem with the default initial
guess of the solvers (all zeros). The results are shown in the column of Table
I named "MPCC_default". As we can see, despite the relaxation of complementary
constraints, the rate of solving the problem successfully with zero warm-start
is almost zero.

Warm Start with Manually Designed Heuristics Another baseline is to warm-start
the formulation with human-designed heuristics. The bookshelf problem requires
inserting one book minimizing the movements of the existing books. Therefore,
a good initial guess is to directly use the continuous and discrete variables
from the original scene. We also pre-solve the separating planes using the
original scene to warm-start the parameters a and $b$. The authors believe
this is the best initial guess readily available without paying much more
effort. The results are shown in Table I in the column "MPCC_Manual".

Implementing this initial guess increases the success rate to more than
$70\%$, significantly better than the default initial guess from the solver.
However, we noticed that this success rate drops as the number of books
increases. When there are only 2 books on the shelf, inserting the 3rd book
has over $80\%$ success rate with the same approach for the initial guess.
Therefore, the manually selected initial guess may not scale to problems with
much larger sizes.

Warm Start with pre-solved Dataset We show the result of using K-nearest
neighbor learner to pick out the top 3 candidates from a 15000 pre-solved
dataset in the column "MPCC_KNN3" column in Table I. The results show that a
simple KNN solver is able to significantly increase the rate of getting a
feasible solution to more than $98\%$. The average iteration of solving the
problem is only slightly more than 1 showing that most problem instances are
solved with only 1 K-nearest neighbor sample. Since the MPCC formulation keeps
the non-convexity of the original formulation, it permits the solver to
explore larger regions, hence the learner does not need to pick out a point or
region that is very close to the optimal solution as required by the convex
formulation. Note that we found warm-starting the non-convex continuous
variables is as important as warm-starting the discrete (complementary)
variables.

###  V-C Mixed-integer Programming Formulation

####  V-C1 Reduced Scale MIP Formulation

We implemented the clustering method proposed in section V-C to reduce the
size of MIP formulation. The dataset used for clustering has a size of 1000,
and is classified into 26 clusters. A random forest classifier is trained to
predict the cluster label on the dataset which achieves 92 %. As the new
problem instance arrives, the random forest classifier selects the cluster,
then the region occupied by the non-convex solutions for this dataset is used
to solve the mixed-integer programming. The results are shown in the columns
named "MIP" in Table I.

As the results show, this method has the highest chance to solve the new
problem instances. However, the solving time is significantly longer than the
MPCC formulation and the convex formulation. The reason is that some problems
which do not get classified correctly. It can take minutes to solve those
cases making the average solving time much longer. We are still working on
tuning the clustering method to eliminate those cases.

####  V-C2 Convex Formulation

To make the online formulation convex, the non-convex part of the continuous
solutions is first discretized, then the complete list of integer variables is
used for training. We discretize the variables using the range given by Table
II. A K-nearest neighbor learner is trained to pick out the first 10 integer
solutions. They are used to turn the mixed-integer program into convex
programs which are then solved one by one until a feasible solution is
retrieved. We select the solver OSQP to solve those convex problems. This
solver automatically uses the previous solution to warm-start the next
solution for faster solving speed. The maximal iterations can be tuned. As
most of the problems take no more than 200 iterations to solve, we set this
parameter to 200. Setting the maximal iteration to 600 will incorporate more
solutions, but the average solving time is longer as it takes more time before
the solver decides infeasibility. The result is shown in the columns named
"Convex" in Table I.

The result shows that the convex formulation achieves the fastest solving
speed due to the strong performance of convex solvers. The authors suspect
that the solvers such as OSQP utilized the heuristics specific to convexity,
hence run faster than well-initialized NLP solvers such as Knitro. On the
other hand, the price to pay for turning the problem into convex is that the
learner needs to perform better than in the MPCC case. This is indicated by
the 4.85 average iterations required to find a feasible solution. Due to the
inherent issue with KNN, some instances are classified incorrectly, and the
solver wastes much longer time going over the infeasible integer variables.

TABLE II: Segmentations of Non-convex Variables variable | range | # of segmentations
---|---|---
Item orientation $\theta$ (rad) | [-pi/2, pi/2] | 8
Separating plane normal $a$ | [-1, 1] | 8
Vertex x position $v_{x}$ (cm) | [-9, 9] | 4
Vertex y position $v_{y}$ (cm) | [0, 11] | 4

As the space is discretized, the $\Theta$ space is no longer Euclidean. The
quantities that are close to each other in the Euclidean distance sense, but
on different sides of the discretization boundary will be assigned different
integer variables, hence no longer close to each other. The KNN utilizes
Euclidean distance, thus ignoring this non-Euclidean effect and failing to
select the correct warm-start integers in those cases.

###  V-D Nonlinear ADMM Formulation

As a non-data-driven benchmark, we implemented the nonlinear ADMM method on
the MINLP formulation by iterating between a mixed-integer formulation and a
nonlinear formulation as described in Section V-D. To keep this method non-
data-driven, manually designed heuristics same as V-B are used to warm-start
the nonlinear formulation to increase the feasibility rate. A sufficient
effort is devoted to tuning the weights G, $\gamma$, and scaling the variables
properly to ensure the performance is on the better end. Since the mixed-
integer approximations on the bilinear constraints induce numerical errors,
termination conditions and accuracy of bilinear collision avoidance
constraints are set lower to match the accuracy of mixed-integer envelope
formulation. The results are shown in the column named "ADMM" in Table I.

According to the results, ADMM works fairly well for getting a feasible
solution. The solving speed is several times slower compared to the data-
driven MPCC formulation and convex formulation. One important reason is that
it usually takes 4-5 iterations (1 iteration includes 1 MIP and 1 NLP) to get
a feasible solution. The nonlinear formulation takes more than 70% of the
solving time, due to the separating plane constraints. If data-driven methods
can be used to warm-start the nonlinear formulation, a faster solving speed
may be achieved. In addition, the average optimal cost from ADMM is worse than
most of the data-driven methods. The authors suspect that the reason is due to
the extra consensus terms in the objective function which guides the
convergence.

###  V-E Hardware Experiment

We implement our optimization results on hardware with a 6 degree of freedom
manipulator [54] to insert a book onto a shelf. To automate the system, a
trajectory planner is required to generate the insertion motion, which is
common in practice. As the focus of this paper is on high level planning, we
simplify this part by manually selecting waypoints. The hardware
implementation is shown in the attached video.

##  VI Conclusion, Discussion and Future Work

This paper proposes data-driven methods to solve MPCC, MIP, or convex
formulation applied to a bookshelf organization problem formulated into mixed-
integer non-convex optimization incorporating mode switch and separating
planes as collision avoidance constraints. The performances of the proposed
data-driven methods are compared. The planned organizations are implemented on
the hardware.

According to our test, data-driven methods can simultaneously achieve good
optimality, fast solving speed, and high success rate leveraging on a
reasonable amount of pre-solved data. Among the online formulations, MPCC
keeps the original non-convex formulations hence can converge better even from
a relatively bad initial warm-start. The implementation of MPCC on NLP solvers
is easier than the convex formulation which requires converting the non-convex
constraints into mixed-integer envelope constraints. The solving speed of MPCC
formulation is also fast. The convex formulation has the fastest solving speed
due to well-performed convex solvers. However, since additional constraints
are added to make the problem convex, it is more difficult for the learner to
select a good warm-start, resulting in more trials. Finally, non-data-driven
methods such as nonlinear ADMM have a decent chance to get a feasible
solution. However, both the optimality and the solving speed are much worse.

The obvious merit of the data-driven approach is that it can solve the problem
fast, but also optimally and reliably as long as the new problem instance is
within the trained parameter regime. An obvious challenge with this method is
how to generalize the training to cases outside the dataset. To achieve this
goal, simple learners such as KNN do not perform well, and more modern
learning methods such as auto-encoders are to be explored.

Acknowledgements Special thanks to Abhishek Cauligi and Professor Bartolomeo
Stellato for helpful discussions, and Yusuke Tanaka and Hyunwoo Nam for
hardware help.

## References

  * [1] R. Deits and R. Tedrake, “Footstep planning on uneven terrain with mixed-integer convex optimization,” in _2014 IEEE-RAS international conference on humanoid robots_. IEEE, 2014, pp. 279–286.
  * [2] X. Lin, J. Zhang, J. Shen, G. Fernandez, and D. W. Hong, “Optimization based motion planning for multi-limbed vertical climbing robots,” in _2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2019, pp. 1918–1925.
  * [3] H. Dai, A. Valenzuela, and R. Tedrake, “Whole-body motion planning with centroidal dynamics and full kinematics,” in _2014 IEEE-RAS International Conference on Humanoid Robots_. IEEE, 2014, pp. 295–302.
  * [4] A. W. Winkler, C. D. Bellicoso, M. Hutter, and J. Buchli, “Gait and trajectory optimization for legged systems through phase-based end-effector parameterization,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 3, pp. 1560–1567, 2018.
  * [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 Robotics and Automation Letters_ , vol. 5, no. 4, pp. 4994–5001, 2020.
  * [6] M. Soler, A. Olivares, P. Bonami, and E. Staffetti, “En-route optimal flight planning constrained to pass through waypoints using minlp,” 2011.
  * [7] S. Boyd and J. Mattingley, “Branch and bound methods,” _Notes for EE364b, Stanford University_ , pp. 2006–07, 2007.
  * [8] J. Tordesillas, B. T. Lopez, and J. P. How, “Faster: Fast and safe trajectory planner for flights in unknown environments,” in _2019 IEEE/RSJ international conference on intelligent robots and systems (IROS)_. IEEE, 2019, pp. 1934–1940.
  * [9] X. Lin, M. S. Ahn, and D. Hong, “Designing multi-stage coupled convex programming with data-driven mccormick envelope relaxations for motion planning,” in _2021 IEEE/RSJ International Conference on Robotics and Automation (ICRA)_. IEEE, 2021.
  * [10] V. Nair, S. Bartunov, F. Gimeno, I. von Glehn, P. Lichocki, I. Lobov, B. O’Donoghue, N. Sonnerat, C. Tjandraatmadja, P. Wang _et al._ , “Solving mixed integer programs using neural networks,” _arXiv preprint arXiv:2012.13349_ , 2020.
  * [11] Y. Tang, S. Agrawal, and Y. Faenza, “Reinforcement learning for integer programming: Learning to cut,” in _International Conference on Machine Learning_. PMLR, 2020, pp. 9367–9376.
  * [12] J.-J. Zhu and G. Martius, “Fast non-parametric learning to accelerate mixed-integer programming for online hybrid model predictive control,” _arXiv preprint arXiv:1911.09214_ , 2019.
  * [13] A. Cauligi, P. Culbertson, E. Schmerling, M. Schwager, B. Stellato, and M. Pavone, “Coco: Online mixed-integer control via supervised learning,” _arXiv preprint arXiv:2107.08143_ , 2021.
  * [14] R. Deits and R. Tedrake, “Computing large convex regions of obstacle-free space through semidefinite programming,” in _Algorithmic foundations of robotics XI_. Springer, 2015, pp. 109–124.
  * [15] L. Dai, Q. Cao, Y. Xia, and Y. Gao, “Distributed mpc for formation of multi-agent systems with collision avoidance and obstacle avoidance,” _Journal of the Franklin Institute_ , vol. 354, no. 4, pp. 2068–2085, 2017\.
  * [16] A. U. Raghunathan, D. K. Jha, and D. Romeres, “Pyrobocop: Python-based robotic control & optimization package for manipulation,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 985–991.
  * [17] M. Stolle and C. G. Atkeson, “Policies based on trajectory libraries,” in _Proceedings 2006 IEEE International Conference on Robotics and Automation, 2006. ICRA 2006._ IEEE, 2006, pp. 3344–3349.
  * [18] C. Liu, C. G. Atkeson, and J. Su, “Biped walking control using a trajectory library,” _Robotica_ , vol. 31, no. 2, pp. 311–322, 2013.
  * [19] Y. Tassa, T. Erez, and W. Smart, “Receding horizon differential dynamic programming,” _Advances in neural information processing systems_ , vol. 20, 2007.
  * [20] A. V. Fiacco, _Introduction to sensitivity and stability analysis in nonlinear programming_. Academic press, 1983, vol. 165.
  * [21] E. N. Pistikopoulos, V. Dua, N. A. Bozinis, A. Bemporad, and M. Morari, “On-line optimization via off-line parametric optimization tools,” _Computers & Chemical Engineering_, vol. 26, no. 2, pp. 175–185, 2002.
  * [22] T. Gal, _Postoptimal Analyses, Parametric Programming, and Related Topics: degeneracy, multicriteria decision making, redundancy_. Walter de Gruyter, 2010.
  * [23] A. Bemporad, M. Morari, V. Dua, and E. N. Pistikopoulos, “The explicit linear quadratic regulator for constrained systems,” _Automatica_ , vol. 38, no. 1, pp. 3–20, 2002.
  * [24] V. Dua and E. N. Pistikopoulos, “Algorithms for the solution of multiparametric mixed-integer nonlinear optimization problems,” _Industrial & engineering chemistry research_, vol. 38, no. 10, pp. 3976–3987, 1999.
  * [25] P. Tøndel, T. A. Johansen, and A. Bemporad, “An algorithm for multi-parametric quadratic programming and explicit mpc solutions,” _Automatica_ , vol. 39, no. 3, pp. 489–497, 2003.
  * [26] V. Dua and E. N. Pistikopoulos, “An algorithm for the solution of multiparametric mixed integer linear programming problems,” _Annals of operations research_ , vol. 99, no. 1, pp. 123–139, 2000.
  * [27] K. Hauser, “Learning the problem-optimum map: Analysis and application to global optimization in robotics,” _IEEE Transactions on Robotics_ , vol. 33, no. 1, pp. 141–152, 2016.
  * [28] G. Tang and K. Hauser, “A data-driven indirect method for nonlinear optimal control,” _Astrodynamics_ , vol. 3, no. 4, pp. 345–359, 2019.
  * [29] T. Power and D. Berenson, “Variational inference mpc using normalizing flows and out-of-distribution projection,” _arXiv preprint arXiv:2205.04667_ , 2022\.
  * [30] H. Dai, G. Izatt, and R. Tedrake, “Global inverse kinematics via mixed-integer convex optimization,” _The International Journal of Robotics Research_ , vol. 38, no. 12-13, pp. 1420–1441, 2019.
  * [31] H. Marchand, A. Martin, R. Weismantel, and L. Wolsey, “Cutting planes in integer and mixed integer,” _Discrete Optimization: The State of the Art_ , vol. 11, p. 397, 2003.
  * [32] H. P. Williams, _Model building in mathematical programming_. John Wiley & Sons, 2013.
  * [33] 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.
  * [34] E. Drumwright, “Rapidly computable viscous friction and no-slip rigid contact models,” _arXiv preprint arXiv:1504.00719_ , 2015.
  * [35] M. Kojima, N. Megiddo, T. Noma, and A. Yoshise, “A unified approach to interior point algorithms for linear complementarity problems: A summary,” _Operations Research Letters_ , vol. 10, no. 5, pp. 247–254, 1991.
  * [36] M. Posa, C. Cantu, and R. Tedrake, “A direct method for trajectory optimization of rigid bodies through contact,” _The International Journal of Robotics Research_ , vol. 33, no. 1, pp. 69–81, 2014.
  * [37] J. Zhang, X. Lin, and D. W. Hong, “Transition motion planning for multi-limbed vertical climbing robots using complementarity constraints,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 2033–2039.
  * [38] S. L. Cleac’h, T. Howell, M. Schwager, and Z. Manchester, “Fast contact-implicit model-predictive control,” _arXiv preprint arXiv:2107.05616_ , 2021.
  * [39] J. Park and S. Boyd, “A semidefinite programming method for integer convex quadratic minimization,” _Optimization Letters_ , vol. 12, no. 3, pp. 499–518, 2018.
  * [40] S. Boyd, S. P. Boyd, and L. Vandenberghe, _Convex optimization_. Cambridge university press, 2004.
  * [41] J. P. Vielma, “Mixed integer linear programming formulation techniques,” _Siam Review_ , vol. 57, no. 1, pp. 3–57, 2015.
  * [42] M. Posa, M. Tobenkin, and R. Tedrake, “Stability analysis and control of rigid-body systems with impacts and friction,” _IEEE Transactions on Automatic Control_ , vol. 61, no. 6, pp. 1423–1437, 2015.
  * [43] P. Bonami and J. Lee, “Bonmin user’s manual,” _Numer Math_ , vol. 4, pp. 1–32, 2007.
  * [44] B. Stellato, G. Banjac, P. Goulart, A. Bemporad, and S. Boyd, “Osqp: An operator splitting solver for quadratic programs,” _Mathematical Programming Computation_ , vol. 12, no. 4, pp. 637–672, 2020.
  * [45] Y. Shirai _et al._ , “Simultaneous contact-rich grasping and locomotion via distributed optimization enabling free-climbing for multi-limbed robots,” _IEEE Proc. 2022 IEEE/RSJ Int. Conf. Intell. Rob. Syst_ , 2022.
  * [46] O. Stein, J. Oldenburg, and W. Marquardt, “Continuous reformulations of discrete–continuous optimization problems,” _Computers & chemical engineering_, vol. 28, no. 10, pp. 1951–1966, 2004.
  * [47] P. M. Castro, “Tightening piecewise mccormick relaxations for bilinear problems,” _Computers & Chemical Engineering_, vol. 72, pp. 300–311, 2015\.
  * [48] X. Lin, G. I. Fernandez, and D. W. Hong, “Reduce: Reformulation of mixed integer programs using data from unsupervised clusters for learning efficient strategies,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 4459–4465.
  * [49] M. Ester, H.-P. Kriegel, J. Sander, and X. Xu, “Density-based spatial clustering of applications with noise,” in _Int. Conf. Knowledge Discovery and Data Mining_ , vol. 240, 1996, p. 6.
  * [50] T. Berthold, “Rens,” _Mathematical Programming Computation_ , vol. 6, no. 1, pp. 33–54, 2014.
  * [51] D. Bertsimas and B. Stellato, “The voice of optimization,” _Machine Learning_ , vol. 110, no. 2, pp. 249–277, 2021.
  * [52] ——, “Online mixed-integer optimization in milliseconds,” _arXiv preprint arXiv:1907.02206_ , 2019.
  * [53] A. Aydinoglu and M. Posa, “Real-time multi-contact model predictive control via admm,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 3414–3421.
  * [54] D. Noh, Y. Liu, F. Rafeedi, H. Nam, K. Gillespie, J.-s. Yi, T. Zhu, Q. Xu, and D. Hong, “Minimal degree of freedom dual-arm manipulation platform with coupling body joint for diverse cooking tasks,” in _2020 17th International Conference on Ubiquitous Robots (UR)_. IEEE, 2020, pp. 225–232.
  * [55] P. Belotti, L. Liberti, A. Lodi, G. Nannicini, A. Tramontani _et al._ , “Disjunctive inequalities: applications and extensions,” _Wiley Encyclopedia of Operations Research and Management Science_ , vol. 2, pp. 1441–1450, 2011.
  * [56] J. P. Vielma and G. L. Nemhauser, “Modeling disjunctive constraints with a logarithmic number of binary variables and constraints,” _Mathematical Programming_ , vol. 128, no. 1, pp. 49–72, 2011.
  * [57] D. Avis and K. Fukuda, “A pivoting algorithm for convex hulls and vertex enumeration of arrangements and polyhedra,” _Discrete & Computational Geometry_, vol. 8, no. 3, pp. 295–313, 1992.

We explain the formulation that is used to solve the MIPs in Sec. V-A.

For a non-convex constraint, we segment it into multiple regions and locally
approximate or relax them into convex constraints. In this paper, we relax the
bilinear constraints locally into convex polytopes (McCormick envelopes). Each
polytope is associated with a unique combination of integer variable values,
hence, mixed-integer convex constraints. Assume the number of regions used is
$N$. Depending on the number of integer variables used, the formulation can
generally be divided into 2 categories: 1) If the number of integer variables
is $N$, we call it $N$ formulation, e.g., the convex hull formulation [55]. 2)
If the number of integer variables is $log_{2}N$, we call it $log_{2}N$
formulation, e.g., [56]. [56] presents a $log_{2}N$ formulation to model the
special ordered sets of type 2 (sos2). However, there are several limitations
of this formulation. For example, the segmented regions need to be connected
to be a valid sos2 constraint, i.e., the two consecutive couple of non-zero
entries can be any consecutive couple in the set. In this paper, we use MIP to
model convex polytopes of arbitrary locations. This can increase the
complexity for sos2 techniques as the polytopes can be disjunctive. The
disjunctive constraints can be handled with convex hull formulations at a
price of introducing more integer variables which may result in slower solving
speeds.

In this appendix, we demonstrate an intuitive but general $log_{2}N$
formulation to model combinations of convex polytopes at any locations which
serves as the base MIP formulation for the bilinear constraints in our paper.
Assume the variable x is enforced to be within one of the $N$ convex
polytopes, denoted by $\textbf{A}_{i}\textbf{x}\leq\textbf{b}_{i}$,
$i=1,...,N$. We introduce $m=log_{2}N$ binary variables $z_{1},...,z_{m}$,
$z_{i}\in\\{0,1\\}$. Each combination of unique values of binary variables can
be assigned to a convex polytope. Let the assignment be:

$\textbf{z}=\bar{\textbf{z}}_{i}\Rightarrow\textbf{A}_{i}\textbf{x}\leq\textbf{b}_{i}$
(6)

Where $\bar{\textbf{z}}_{i}=[\bar{z}_{i,1},...,\bar{z}_{i,m}]$ are constant
binary values associated with polytope $i$. Note $\bar{z}_{i}\neq\bar{z}_{j}$
if $i\neq j$. In other words, we require that when
$\textbf{z}=\bar{\textbf{z}}_{i}$, x stays within the polytope
$\textbf{A}_{i}\textbf{x}\leq\textbf{b}_{i}$; otherwise, the constraint is
unenforced.

Denote the vertices of the polytopes by
$\textbf{v}_{i,1},...,\textbf{v}_{i,n_{i}}$, $i=1,...,N$, where $n_{i}$ is the
number of vertices associated with polytope $i$. Each vertex
$\textbf{v}_{i,j}$ is assigned a continuous non-negative variable
$\lambda_{i,j}\in[0,1]$. In general, one can run a mathematical program (e.g.
[57]) to get vertices from the nondegenerate system of inequalities
$\textbf{A}_{i}\textbf{x}\leq\textbf{b}_{i}$. As a result, the assignment
becomes:

$\textbf{z}=\bar{\textbf{z}}_{i}\Rightarrow\ \ \begin{aligned}
&\textbf{x}=\sum_{j=1}^{n_{i}}\lambda_{i,j}\textbf{v}_{i,j}\\\
&\sum_{j=1}^{n_{i}}\lambda_{i,j}=1,\ \ \lambda_{i,j}\in[0,1]\end{aligned}$ (7)

The formulation can be written as:

$\displaystyle(5.a)\quad\textbf{x}=\sum_{i=1}^{N}\sum_{j=1}^{n_{i}}\lambda_{i,j}\textbf{v}_{i,j}$
(8) $\displaystyle(5.b)\quad\sum_{i=1}^{N}\sum_{j=1}^{n_{i}}\lambda_{i,j}=1,\
\ \lambda_{i,j}\in[0,1]$ $\displaystyle(5.c)\quad\sum_{k=1,...,N}^{k\neq
i}\sum_{j=1}^{n_{k}}\lambda_{k,j}\leq\sum_{l=1}^{m}|z_{l}-\bar{z}_{i,l}|$

The set of constraint in (8) enforces (7) for all polytopes, as
$\sum_{l=1}^{m}|z_{l}-\bar{z}_{i,l}|=0$ only when
$\textbf{z}=\bar{\textbf{z}}_{i}$, enforcing that all $\lambda$’s that are not
associated with polytope $i$ to be zero. If
$\textbf{z}\neq\bar{\textbf{z}}_{i}$, $\sum_{l=1}^{m}|z_{l}-\bar{z}_{i,l}|\geq
1$ and constraint (5.c) is looser than constraint (5.b), hence, trivial.

Note that formulation (8) works for any convex polytope that can be written as
$\textbf{A}_{i}\textbf{x}\leq\textbf{b}_{i}$. In this paper, the polytopes are
McCormick envelope constraints which is a special case.


††thanks<EMAIL_ADDRESS>
# Non volatile photo-switch using a diamond pn junction

Cédric Masante  Univ. Grenoble Alpes, CNRS, Institut Néel, 38000 Grenoble,
France  Martin Kah  Univ. Grenoble Alpes, CNRS, Institut Néel, 38000 Grenoble,
France  Clément Hébert  Univ. Grenoble Alpes, INSERM, U1213, Grenoble Institut
des Neurosciences (GIN), Grenoble, France  Nicolas Rouger  LAPLACE, Université
de Toulouse, CNRS, INPT, UPS, Toulouse, France  Julien Pernot  Univ. Grenoble
Alpes, CNRS, Institut Néel, 38000 Grenoble, France

###### Abstract

Ultrawide bandgap semiconductor technologies offer potentially revolutionary
advances in the rapidly developing areas of quantum communication, short
wavelength optics, smart energy conversion and biomedical interfaces. These
strongly demanding technologies can be partly constructed using conventional
devices but new hybrid architectures are needed to overpass current
performances and add functionalities. Here, we propose a new concept based on
the specific properties of a diamond pn junction combined with both an
electric and optical control of the depletion region. Using this junction as a
gate in a junction field effect transistor, we report a proof of concept of a
non volatile diamond photo-switch. A diamond pn junction made with nitrogen
deep donors in the n?side is demonstrated to be optically activated thanks to
visible light. The n-type diamond gate is almost devoid of free carriers in
the dark and thus insulating. Illuminating the device renders the standard
electrical gate control of the transistor efficient. Without illumination, the
device is frozen, keeping a permanent memory of the current state. This new
way of operating the device opens numerous possibilities to store and transfer
information or energy with applications in the field of electrical aircraft or
aerospace electronics, power electronics, bio-electronics and quantum
communication.

##  I Introduction

Because of their exceptional characteristics, ultrawide bandgap semiconductors
(UWBGS) like diamond, aluminium nitride, boron nitride, gallium oxide are
already in focus for producing devices beyond state-of-the-art semiconductor
technologies Tsao2018 . For electronics, most researchers and engineers are
still using concepts envisaged for silicon and wide bandgap semiconductors
(like silicon carbide and gallium nitride) where much smaller band gap
energies are at stake. The limits of this incremental approach for the
fabrication of UWBGS-based devices are already well known and established. In
this paper, we propose to go beyond current paradigms by putting forward the
design of novel hybrid technologies (optical and electrical) based on
unconventional concepts which exploit the huge energy scale of UWBGS.

Diamond’s outstanding properties like wide bandgap, elevated carrier mobility,
high breakdown field and huge thermal conductivity make it a unique material
for high power and/or high temperature and/or radiation-resistant electronics.
The long valley-polarization relaxation times of electrons make diamond a good
candidate for use as valleytronics transistor IsbergNatMater2013 ;
Suntornwipat2021 . Diamond takes also advantage of point defects, like
nitrogen vacancy centres, for applications in quantum sensing Maminscience2013
, quantum communication HensenNature2015  or spin-based photonic quantum
technologies AtatureNatureReviews2018 . Finally, its bio-compatibility offers
the possibility of integrating diamond electronics into the human body
Rifai2021 ; Lichter2015 ; Fox2020 ; Mani2020 ; Fan2020 ; Dankerl2009 ;
Piret2015 . To summarize, this material provides new opportunities to locally
probe electric or magnetic signals, to reliably process data and to safely
communicate them in harsh or living environments. Diamond devices capable of
storing charges in a non volatile manner or addressing individual devices are
needed to develop monolithic diamond-based systems. As examples, the charge
accumulations could be used for power supply (local on-chip energy storage) or
information storage with read and write access (non volatile memory) and
individual on-chip addressing for application in neuronal multiplexers.
Embedded on-chip devices coupled with an ultra high conversion efficiency of
beta-voltaic cell using diamond $pn$ junction ShimaokaAPL2020  can be
envisaged for aerospace or harsh environment applications.

In this paper, we report on a diamond-based non volatile photo-switch with a
permanent off state having no loss. The concept is demonstrated by using
monocristalline $p-$type diamond deposited by plasma enhanced chemical vapour
deposition on a highly nitrogen-doped diamond substrate manufactured using
high pressure and high temperature. In the first part, the concept will be
described. Then, the active and freeze state of the $pn$ junction will be
established and the proof of concept of a non volatile photo-switch will be
demonstrated thanks to an electro-optical control of a diamond junction field
effect transistor. Then, we will discuss examples of application fields in
which this concept can efficiently operate such as power electronics, quantum
communication or bio-electronics.

##  II Non volatile diamond-based photo-switch concept

For analogue or digital operations, the ubiquitous silicon-based field effect
transistors work by using a gate bias voltage to deplete (OFF state) or
populate (ON state) the channel. In the hybrid structure proposed here, the
external illumination activates or freezes the device. In this way, the device
can be driven by a combination of electrical and optical signals, providing a
higher degree of freedom. In conventional semiconductors with shallow energy
levels associated to donors or acceptors, doping level and carrier
concentrations are almost equal at room temperature. The doped material is
conducting. In nitrogen doped diamond, nitrogen impurities in substitutional
sites create deep donor levels at 1.7 eV from the conduction band.
Consequently, at room temperature, electrons remain bound to the dopants. Even
at high doping levels (1018 cm-3 to 1019 cm-3) nitrogen doped (N-doped)
diamond samples are completely insulating, and so cannot be electrically
driven. This is the freeze state of the device. Thanks to external
illumination, the absorption of photons promotes some of the bound electrons
to the conduction band, increasing the free carrier concentration. The device
thus becomes conducting and prone to be electrically driven. This is the
active state of the device. The control of the ON and OFF states of the
transistor is still ensured by a gate bias voltage. However, the switch from
ON to OFF will only be achievable during the active state. This concept has to
be distinguished from photo-conductive semiconductor switch devices majda-
zdancewiczCurrentStatePhotoconductive2018  which have also been demonstrated
to use N-doped diamonds hallPhotoconductiveSwitchHigh2020  well. Optically
generated carriers most participate in conduction which means these devices
require a high power pulsed laser source to be efficient. For example, the
device reported in hallPhotoconductiveSwitchHigh2020  makes use of a pulsed
laser with pulse power of about 6.7 kW (the average power is 25 mW). In this
work, a bulk boron doped ($p-$type) diamond is used as the channel and the
nitrogen doped diamond as the gate. The lower ionization energy of boron
acceptor (0.38 eV) means it can conduct at room temperature in the dark while
the gate requires external optical excitation to generate free electrons.
Optically generated carriers do not significantly participate in the channel
conduction. Only the n-type gate junction is optically controlled by the light
source.

##  III Optically controlled diamond pn junction

A schematic description of the sample structure is shown in Figure 1a and top
view images in Figure 1b and 1c. A white LED driven by a tuneable current was
used to optically excite the sample from the top surface. The average power
density of the beam was measured with a photo-diode detector, ranging from 0
to 11 mW.cm-2. Part of the light was absorbed and reflected by the top
metallic layer. The Figure 1d shows a room temperature current-voltage
characteristic under various optical power varying between 0 and 11 mW.cm-2. A
diode characteristic was clearly measured under optical excitation while the
dark current was below detection limit (< 10-8 A.cm-2) in the whole voltage
range. The threshold voltage was around 3 V, as expected for a diamond $pn$
junction where the Fermi level lies close to the deep nitrogen donor level in
n-side region. This clearly demonstrates that a nitrogen doped and boron doped
diamond junction can be optically activated by continuous low power white
light excitation. The LED does not emit in the UV range required for electron-
hole generation in diamond ($\lambda<225$ nm) KoizumiSceince2001 . Also, no
intrinsic carriers were generated from the valence band to the conduction
band. The linear relationship between the substrate conductivity and the
irradiance (Figure 1d) is characteristic of a one photon impurity ionization
process. In this experiment, the white light was found to be efficient to
observe a photo-current because the transition from the nitrogen impurity
level to the conduction band is efficient for wavelengths lower than 550 nm
with absorption coefficients of the order of 10 to 100 cm-1
deweerdtDeterminationDefectConcentration2008 . The corresponding absorption
depths were between 1 and 0.1 mm, from the diamond surface. As expected, the
forward current was greatly limited by the thick, resistive $n-$type substrate
with a current density in the $\mu$A.cm-2 range in the forward regime at room
temperature (the current was normalized by the $p-$type layer surface area,
delineated by the mesa etching). One can however extrapolate that the pn
junction conductivity could be increased by orders of magnitude by (i) using a
higher power continuous light sources up to 1000 W.cm-2 light output power
than can be achieved with LEDs or lasers, (ii) depositing transparent top
contacts or illuminating from the backside of the sample, (iii) using a <500
nm light source to maximize the optical absorption, (iv) using a thinner
N-doped substrate and (vi) forming a proper ohmic contact on the N-doped
diamond. However this simple setup provides a strong proof of concept of the
nitrogen-doped diamond electrical activation under light exposure.

Figure 1: Optically controlled diamond pn junction a Schematic representation
of the device and the measurement setup during the $pn$ diode electrical
characterization. The $p-$type layer resistance is neglected in comparison
with the substrate resistance. b Optical view of the sample. c  Secondary
electron microscope top view image of a singular test device (colours added).
d Room temperature current-voltage characteristic of the drain-gate $pn$ diode
activated by a white light LED with various optical power ranging between 0
and 11 mW/cm?. The source is floating and the gate/substrate is grounded. The
current is normalized by the $p-$type layer area delineated by mesa etching. e
Current-voltage characteristic of the $pn$ diode under light exposure at
$-50\text{\,}\mathrm{\SIUnitSymbolCelsius}$.

##  IV High breakdown capability of the pn junction gate

As an UWBGS, diamond is expected to be able to sustain an extreme electric
field before reaching its breakdown. Hiraiwa et
al.hiraiwaFigureMeritDiamond2013  estimated the diamond breakdown field for a
doping level [NA]$=2\times 10^{17}$ cm-3 at 10 MV.cm-1 and Kamakura et
al.kamakuraInitioStudyAvalanche2016  at 7 MV.cm-1, compared to 0.2 to 0.6
MV.cm-1 for silicon. However experimental demonstrations are challenging in
diamond devices at the moment due to the presence of a strong leakage current
which causes an early soft breakdown. To circumvent this difficulty in this
original structure, we performed the measurement at
$-50\text{\,}\mathrm{\SIUnitSymbolCelsius}$ in order to reduce the temperature
dependent leakage current. Since the substrate carriers were optically
generated, and the substrate was still more resistive than the $p-$layer, no
significant variation of the diode resistance was measured. In the dark, the
reverse bias is held by the full thickness of the insulating substrate and the
current is below the detection limit for a bias up to 1 kV. However when the
diode was reverse biased under optical excitation, because the N-doped
substrate is activated with a doping level much higher than the $p-$type
layer, the Space Charge Region (SCR) was found to extend almost exclusively in
the $p-$side. Almost all the potential drop is then sustained by the 400 nm
boron doped layer. The Figure 1e shows the obtained IV curve, exhibiting a
clear change of slope at V=-156 V attributed to the junction’s breakdown. Due
to the large substrate serial resistance, the breakdown characteristic did not
exhibit the sharp increase of current with a power law, but was instead
immediately limited by the substrate resistance and appears linear. A local
breakdown, not occurring on the full surface of the diode, might explain the
higher serial resistance measured after the breakdown compared to the forward
regime. This is represented by an additional serial resistance $R_{BD}$ in
Figure 1e. The maximum electric field at breakdown has been estimated to be
5.2 MV.cm-1 in a 1D punch-through configuration. This is a remarkably high
value which is among the highest critical fields reported for diamond and much
larger than other common wide band gap materials such as SiC and GaN (ideal
breakdown around 2-3 MV.cm-1). Since no special care to reduce the 2D and 3D
effects was implemented in the fabrication process, the real value is most
probably higher due to the well known electric field crowding effect. The high
breakdown voltage of diamond pn junction will ensure the fabrication of
reliable and robust devices able to work in high voltage operations under
illumination but even more importantly, it ensures a galvanometric insulation
of the source to the gate of the junction field effect transistor described
below without illumination.

##  V Electro-optical control of a diamond transistor

A junction field effect transistor device was fabricated on the same sample to
demonstrate the first non volatile photo-switch using diamond pn junction. The
current flows through the bulk of the $p-$type layer between the drain and
source under bias. Under optical illumination, the gate contact on the back
side of the substrate controls the SCR extension in the channel. The drain to
source distance is 10 $\mu$m and contacts have a 70 $\mu$m width. It is
remarkable that room temperature transistor characteristics under an optical
power of 11 mW/cm? were obtained as shown in Figure 2a and 2b. The typical ON-
state linear regime, pinch-off and OFF-state are schematically drawn in Figure
2c, 2d, and 2e, respectively. The pinch-off is due to the extension of the SCR
close to the drain contact, when the drain to gate bias approaches or is
greater than the threshold bias (i.e. the drain to source bias $V_{DS}$ is
around -30 V or more for $V_{SUB}=0$ V). Due to the low carrier density in the
nitrogen doped substrate, the gate junction can be forward biased (i.e.
$V_{SUB}\lesssim-4$ V) without any significant gate current. The SCR modulated
by the gate bias can be fully removed to increase the channel conduction (0.25
mA/mm). Figure 2a shows that a gate to source bias $V_{SUB}=-10$ V
significantly increases the drain current. The gate current remains orders of
magnitude lower, i.e. less than $2\times 10^{-6}$ mA.mm-1. No drain current
modulation from the gate could be observed without optical illumination.

Figure 2: Non volatile photo-switch using a diamond junction field effect
transistor. a Quasi-static room temperature transistor characteristic, drain
current versus drain-source voltage, measured for a substrate to source bias
VSub varying from -10 V to +30 V with 10 V step and continuous 11 mW/cm?
optical power density. b Drain current density as function of VSub in the
Pinch-off regime at VDS=-40 V. c Schematic cross section of the device in ON-
state at low VDS (linear regime), d at high VDS (pinch-off regime) and e in
OFF-state. f Measured switched drain current over one commutation cycle of the
fabricated JFET at $V_{DS}=-40$ V, as function of the light and substrate
bias. The device is only able to commute under an optical signal.

A measured commutation cycle at a fixed drain to source bias V${}_{DS}=-40V$
(Pinch-off) is shown in Figure 2f, where the freezing of the device’s
electrostatic can be optically controlled. The following paragraph describes
each step represented in Figure 2f: (1), the substrate to source bias VSub is
set at 0 V with the light on, a bulk conduction channel between the drain and
source is obtained (normally-on device). (2), the light is switched off, the
free electron concentration in the N-doped substrate is extremely low (freeze
state). In an ideal case without compensation from acceptor ions at room
temperature, the free electron concentration is below $10^{4}$ $cm^{-3}$. (3)
The gate substrate electric command was applied (V${}_{Sub}=40V$) to fully
deplete the channel (OFF state). However, because the N-doped gate is devoid
of free carriers in the dark, no electrostatic modulation can occur in the
dark and the channel remains conductive even if the gate voltage corresponds
to a OFF state. (4) The light is switched on, providing the optical signal
activating the device (active state). The device was able to respond to the
electrical gate command and switches to OFF state. The channel was fully
depleted and the drain current $I_{DS}$ dropped below the detection limit. (5)
The light was switched off (freeze state). (6) The gate electric command to
switch the device in ON state was applied, V${}_{Sub}=0$ V, but the channel
remained insulating due to the absence of optical excitation. (7) The light
was switched on (active state) and the device was able to switch to the ON
state.

This sequence was replicated several times as shown in Figure 3a to ensure the
reliability of the measurements. Although the first cycle seems different, all
the following cycles were found to exhibit the same characteristic. In order
to check the non volatile character of this optically activated switch, the
device was tested during a long time in the off and freeze state. Figure 3b
and 3c shows the drain current versus time for VDS=-40 V at room temperature
and $100\text{\,}\mathrm{\SIUnitSymbolCelsius}$, respectively for long
duration (2 days at room temperature and 140 minutes for 373 K). One can
observe that the $I_{ON}/I_{OFF}$ ratio is larger than 4 orders of magnitude
during more than 48 hours at room temperature and larger than 3 orders of
magnitude during more than 1 hour at
$100\text{\,}\mathrm{\SIUnitSymbolCelsius}$. This non-volatile photo-switch
benefits from the unique properties of diamond such as its huge breakdown
field and high temperature operation.

The large RC time constant and, more precisely, the huge resistance R caused
by the high substrate’s thickness and resistivity meant that the turn-on time
was long but the modulation was clearly observed. Effects of traps most likely
present at the substrate/epilayer and epilayer/oxide interfaces could not be
distinguished. Reducing the turn-on time is a key factor in obtaining an
efficient device, but as discussed previously, this could be decreased by
orders of magnitude as it directly depends on the N-doped layer resistance.
Additionally, in this kind of FET device the optical excitation only has to be
present during the short time of the device switching. As well as the
previously mentioned ways to increase this simple proof of concept
performances, the use of a pulsed light source could provide a stronger
optical excitation with limited total power usage.

Figure 3: Temporal stability of the non-volatile photo-switch.  a Drain
current density versus time for VDS=-40 V at T=297 K for multiple commutation
cycles alternating freeze/active and on/off states by using various substrate
voltages and optical powers as described by the chronogram. Drain current
versus time for VDS=-40 V b at T= 297 K and c T= 373 K , for long duration in
order to illustrate the non volatile character of the off state in the freeze
state.

##  VI Non volatile photo-switch for new applications

The next challenge to create new devices will be to grow nitride-based light
emitting diodes (LED) directly on diamond transistor. Indeed, creating
electrically driven visible or sub-band gap UV LEDs on top of the diamond
transistors will enable high optical power density excitation to be attained.
This would in turn enable the efficient control of the non volatile photo-
switch, i.e. to read or write the transistor state (ON or OFF thanks to the
FET gate) and store or reset the transistor state (freeze or activate thanks
to the LED electrode). This new device paves the way for numerous practical
applications where the superior properties of diamond semiconductors can be
advantageously exploited. For embedded electronics, non volatile memories,
i.e. the association of a switch with a capacitance for example, could be very
useful to store and read information in a harsh environments (radiation
resistant and high temperature) like in aerospace and aeronautic applications.
For power electronics, the gate terminal of classical power transistors is
affected by high $dV/dt$ and Miller coupling, causing undesired turn-ON and
possible catastrophic shoot-through. Thanks to the perfect galvanic insulation
of the gate terminal without illuminating light, the Miller coupling between
the Drain and the Gate is completely removed thus offering the strongest
immunity to parasitic turn ON. For quantum communication, the LED will make
the excitation of a nitrogen vacancy (NV) centre possible but more importantly
the diamond FET will also permit control of electrostatics and a non-volatile
storage of the NV charge state. Finally, we shall describe a bio-electronics
applications in which the unique properties of diamond in the context of this
concept are fully exploited. Diamond-based electronic is a perfect candidate
for chronic biomedical implants since it provides solution for most of the
scientific barriers that prevent most advanced technology to be transferred to
clinical applications. Its chemical inertness combined with its carbon nature
make it a biocompatible material compatible with all the standard medical
imaging techniques including Magneto Resonance Imaging Hebert2015 .
Furthermore, diamond electronic can be easily sterilized without damaging its
performances. Diamond electrodes and transistors have thus already proved to
be very interesting candidate combination for interaction with the body for
both sensing and stimulation Dankerl2009 ; Piret2015 . This is particularly
true in the case of neural implants where bidirectional communication is
required to develop neural prosthesis, help neural regeneration and gain
further understanding of brain signal processing and dynamics. Future diamond
embedded electronics have the potential to create implants with a high level
of stability that do not require heavy hermetic insulation as is the case with
silicon and other emerging polymer and flexible logic electronics Mani2020 .
The concept offered by the photo-switch using a diamond pn junction is one of
the first diamond implantable logic electronics which can be combined with a
diamond sensor to develop a low power diamond time division multiplexer. Each
line of an array of solution gated field effect diamond transistor can be
sequentially read as in the case of CCD sensor, allowing for the reduction of
the cumbersome electric wires that dramatically increase the size of the
implant Schaefer2020 . Additional the freezing state of the switch means low
power consumption can be achieved which is a key parameter for chronic
implants.

##  VII Conclusion

We demonstrated a novel concept of device with a hybrid optical and electric
gate control. The $pn$ junction made of N-doped diamond can be optically
controlled using a white light LED. A breakdown of -156 V at
$-50\text{\,}\mathrm{\SIUnitSymbolCelsius}$ has been measured, corresponding
to a peak electric field estimated to be superior to 5.2 MV.cm-1. This proof
of concept of hybrid electrical and optical switching device consists of a
diamond $p-$type channel JFET. A clear transistor effect was measured under
light exposure (active state) while almost no current modulation was observed
in the dark (freeze state). The robustness of the non volatile photo-switch
has been demonstrated during two days at room temperature and more than 140
minutes at $100\text{\,}\mathrm{\SIUnitSymbolCelsius}$ The high turn-on and
turn-off times are the consequence of the thick, resistive, substrate used as
the gate, the low light power density and the opaque metallic top contacts.
They could be drastically reduced by simple changes to design of the device,
opening the route to numerous applications in the field of power electronics,
bio-electronics and quantum devices.

##  VIII Experimental Section

###  VIII.1 Diamond pn junction structure deposition

The sample structure is depicted in Figure 1a. The $p-$type layer consists of
a 400 nm boron-doped ([NA-ND]$=2\times 10^{17}$ cm-3) diamond layer grown by
microwave plasma assisted CVD, provided by DiamFab, on a $4\times 4\times 0.5$
mm3 Sumitomo High Pressure High Temperature (HPHT) Nitrogen-doped ($n-$type)
diamond substrate, also referred as 1b diamond. Its nitrogen content is
estimated to be in the order of $10^{19}$ cm-3, mostly in single
substitutional sites. The precise doping level was not measured but the sample
exhibited the characteristic yellow tint of a high nitrogen content, as shown
in the optical top view of the sample in Figure 1b. But, conversely, the
$p-$type layer doping level was confirmed by capacitance-voltage measurements,
not shown here, performed on metal oxide semiconductor capacitors test
devices, fabricated on the same sample.

###  VIII.2 Fabrication of diamond transistor

To lower the energy barrier between the top metal contacts and the $p-$type
epilayer, a 200 nm highly boron-doped ([B]$>4\times 10^{20}$ cm-3) diamond
layer was selectively grown using an aluminium mask. A mesa etching was then
performed down to the substrate by reactive ion etching with an oxygen plasma.
This enabled the electric separation of each test devices on the sample
surface from each other. The top metallic contacts were formed by a
Titanium/Platinum/Gold (30 nm/20 nm/30 nm) stack deposited by electron beam
evaporation and annealed for 30 min at
$500\text{\,}\mathrm{\SIUnitSymbolCelsius}$ in high vacuum condition
($<10^{-7}$ mBar). An UV ozone treatment was performed on the surface before
the deposition of 50 nm alumina layer by atomic layer deposition to passivate
the whole surface, followed by post annealing in high vacuum condition
($<10^{-7}$ mBar) at $500\text{\,}\mathrm{\SIUnitSymbolCelsius}$ for 1 hour
Pham2018 . The top drain and source contact consists of a
Titanium/Platinum/Gold (30 nm/20 nm/30 nm) stack, deposited by electron beam
evaporation on the 50 nm alumina layer. No lithography step has been performed
after the passivation oxide deposition. Also the metal oxide metal stack thus
formed is electrically broken at the time of the first measurement by applying
a few volts bias voltage. The detailed process is given in Ref. Masante2021 .
The back contact between the substrate and the gold layer on the alumina plate
was done thanks to silver paste.

###  VIII.3 Electrical characterization of the device

The PN diode electrical characterization was performed using a Keithley 2612.
The JFET electrical characterization was performed using a Modulab Solartron
including the low current module. The electrical contacts were obtained using
probes moved by a piezoelectric system with a LINKAM system to control
temperature under vacuum. The microscope white LED light beam was used to
optically excite the sample from the top. The use of a nitrogen-doped
substrate as $n-$type layer is convenient and greatly simplifies the
fabrication process. However it also causes a very large serial resistance
$R_{Sub}$ due to its thickness (0.5 mm) and large resistivity. Also the silver
paste back contact is probably of poor quality. Capacitance measurements are
prohibited by the low cut-off frequency induced by the serial resistance so we
only relied on static electrical measurements. This simple design was however
found to be well suited to demonstrate this proof of concept.

###### Acknowledgements.

The authors are thankful to David Eon (Univ. Grenoble Alpes, CNRS, Institut
Néel, France), Etienne Gheeraert (Univ. Grenoble Alpes, CNRS, Institut Néel,
France), Juliette Letellier (DiamFab, France), Khaled Driche (DiamFab,
France), Hitoshi Umezawa (National Institute of Advanced Industrial Science
and Technology, Japan) and Robert Nemanich (Arizona State University, USA) for
fruitful discussions.

## References

  * (1) J. Y. Tsao, S. Chowdhury, M. A. Hollis, D. Jena, N. M. Johnson, K. A. Jones, R. J. Kaplar, S. Rajan, C. G. V. de Walle, E. Bellotti, C. L. Chua, R. Collazo, M. E. Coltrin, J. A. Cooper, K. R. Evans, S. Graham, T. A. Grotjohn, E. R. Heller, M. Higashiwaki, M. S. Islam, P. W. Juodawlkis, M. A. Khan, A. D. Koehler, J. H. Leach, U. K. Mishra, R. J. Nemanich, R. C. N. Pilawa?Podgurski, J. B. Shealy, Z. Sitar, M. J. Tadjer, A. F. Witulski, M. Wraback, J. A. Simmons, Advanced Electronic Materials 2018, 4, 1600501.
  * (2) J. Isberg, M. Gabrysch, J. Hammersberg, S. Majdi, K. K. Kovi, D. J. Twitchen, Nat Mater 2013, 12, 760.
  * (3) N. Suntornwipat, S. Majdi, M. Gabrysch, K. K. Kovi, V. Djurberg, I. Friel, D. J. Twitchen, J. Isberg, Nano Lett. 2021, 21, 868.
  * (4) H. J. Mamin, M. Kim, M. H. Sherwood, C. T. Rettner, K. Ohno, D. D. Awschalom, D. Rugar, Science 2013, 339, 557.
  * (5) B. Hensen, H. Bernien, A. E. Dréau, A. Reiserer, N. Kalb, M. S. Blok, J. Ruitenberg, R. F. L. Vermeulen, R. N. Schouten, C. Abellán, W. Amaya, V. Pruneri, M. W. Mitchell, M. Markham, D. J. Twitchen, D. Elkouss, S. Wehner, T. H. Taminiau, R. Hanson, Nature 2015, 526, 682.
  * (6) M. Atature, D. Englund, N. Vamivakas, S.-Y. Lee, J. Wrachtrup, Nature Reviews Materials 2018, 3, 38.
  * (7) A. Rifai, S. Houshyar, K. Fox, Annals of 3D Printed Medicine 2021, 1, 100002.
  * (8) S. G. Lichter, M. C. Escudié, A. D. Stacey, K. Ganesan, K. Fox, A. Ahnood, N. V. Apollo, D. C. Kua, A. Z. Lee, C. McGowan, A. L. Saunders, O. Burns, D. A. X. Nayagam, R. A. Williams, D. J. Garrett, H. Meffin, S. Prawer, Biomaterials 2015, 53, 464.
  * (9) K. Fox, N. Mani, A. Rifai, P. Reineck, A. Jones, P. A. Tran, A. Ramezannejad, M. Brandt, B. C. Gibson, A. D. Greentree, N. Tran, ACS Appl. Bio Mater. 2020, 3, 29.
  * (10) N. Mani, A. Rifai, S. Houshyar, M. A. Booth, K. Fox, Medical Devices and Sensors 2020, 3, e10127.
  * (11) B. Fan, C. A. Rusinek, C. H. Thompson, M. Setien, Y. Guo, R. Rechenberg, Y. Gong, A. J. Weber, M. F. Becker, E. Purcell, W. Li, Microsystems and Nanoengineering 2020, 6, 1.
  * (12) M. Dankerl, S. Eick, B. Hofmann, M. Hauf, S. Ingebrandt, A. Offenhausser, M. Stutzmann, J. A. Garrido, Advanced Functional Materials 2009, 19, 2915.
  * (13) G. Piret, C. Hébert, J.-P. Mazellier, L. Rousseau, E. Scorsone, M. Cottance, G. Lissorgues, M. O. Heuschkel, S. Picaud, P. Bergonzo, B. Yvert, Biomaterials 2015, 53, 173.
  * (14) P. Siyushev, M. Nesladek, E. Bourgeois, M. Gulka, J. Hruby, T. Yamamoto, M. Trupke, T. Teraji, J. Isoya, F. Jelezko, Science 2019, 363, 728.
  * (15) T. Shimaoka, H. Umezawa, K. Ichikawa, J. Pernot, S. Koizumi, Appl. Phys. Lett. 2020, 117, 103902.
  * (16) E. Majda-Zdancewicz, M. Suproniuk, M. Paw?owski, M. Wierzbowski, Opto-Electronics Review 2018, 26, 92.
  * (17) D. L. Hall, L. F. Voss, P. Grivickas, M. Bora, A. M. Conway, P. ??ajev, V. Grivickas, IEEE Electron Device Letters 2020, 41, 1070.
  * (18) S. Koizumi, K. Watanabe, M. Hasegawa, H. Kanda, Science 2001, 292, 1899.
  * (19) F. De Weerdt, A. T. Collins, Diamond and Related Materials 2008, 17, 171.
  * (20) A. Hiraiwa, H. Kawarada, Journal of Applied Physics 2013, 114, 034506.
  * (21) Y. Kamakura, T. Kotani, K. Konaga, N. Minamitani, G. Wakimura, N. Mori, in 2015 IEEE International Electron Devices Meeting (IEDM), 2015, p. 5.2.1-5.2.4.
  * (22) C. Hébert, J. Warnking, A. Depaulis, L. A. Garcon, M. Mermoux, D. Eon, P. Mailley, F. Omnès, Materials Science and Engineering: C 2015, 46, 25.
  * (23) N. Schaefer, R. Garcia-Cortadella, J. Martinez-Aguilar, G. Schwesig, X. Illa, A. M. Lara, S. Santiago, C. Hébert, G. Guirado, R. Villa, A. Sirota, A. Guimera-Brunet, J. A. Garrido, 2D Mater. 2020, 7, 025046.
  * (24) T. T. Pham, M. Gutiérrez, C. Masante, N. Rouger, D. Eon, E. Gheeraert, D. Araujo, J. Pernot, Appl. Phys. Lett. 2018, 112, 102103.
  * (25) C. Masante, N. Rouger, J. Pernot, J. Phys. D: Appl. Phys. 2021, 54, 233002.