# Updated Core Libraries of the ALPS Project

Markus Wallerberger<EMAIL_ADDRESS>Sergei Iskakov  Alexander Gaenko
Joseph Kleinhenz  Igor Krivenko  Ryan Levy  Jia Li  Hiroshi Shinaoka  Synge
Todo  Tianran Chen  Xi Chen  James P. F. LeBlanc  Joseph E. Paki  Hanna
Terletska  Matthias Troyer  Emanuel Gull  University of Michigan, Ann Arbor,
Michigan 48109, USA  University of Illinois, Urbana-Champaign, Illinois 61820,
USA  Department of Physics, Saitama University, Saitama 338-8570, Japan
Department of Physics, University of Tokyo, Tokyo 113-0033, Japan  Institute
for Solid State Physics, University of Tokyo, Kashiwa 277-8581, Japan  West
Chester University of Pennsylvania, West Chester, Pennsylvania 19383, USA
Center for Computational Quantum Physics, Flatiron Institute, New York, NY
10010, USA  Department of Physics and Physical Oceanography, Memorial
University of Newfoundland, St. John’s, Newfoundland & Labrador A1B 3X7,
Canada  Middle Tennessee State University, Murfreesboro, Tennessee 37132, USA
Institut für Theoretische Physik, ETH Zurich, 8093 Zurich, Switzerland

###### Abstract

The open source ALPS (Algorithms and Libraries for Physics Simulations)
project provides a collection of physics libraries and applications, with a
focus on simulations of lattice models and strongly correlated electron
systems. The libraries provide a convenient set of well-documented and
reusable components for developing condensed matter physics simulation codes,
and the applications strive to make commonly used and proven computational
algorithms available to a non-expert community. In this paper we present an
update of the core ALPS libraries. We present in particular new Monte Carlo
libraries and new Green’s function libraries.

††journal: Computer Physics Communications

PROGRAM SUMMARY

Program Title: ALPS Core libraries 2.3
Project homepage: http://alpscore.org
Catalogue identifier: –
Journal Reference: –
Operating system: Unix, Linux, macOS
Programming language: `C++11` or later
Computers: any architecture with suitable compilers including PCs, clusters
and supercomputers
Licensing provisions: GNU General Public License (GPLv2+)
Classification: 6.5, 7.3, 20
Dependencies: CMake ($\geq 3.1$), Boost ($\geq 1.56$), HDF5 ($\geq 1.8$),
Eigen ($\geq 3.3.4$)
Optional dependencies: MPI ($\geq 2.1$), Doxygen
Nature of problem: Modern, lightweight, tested and documented libraries
covering the basic requirements of rapid development of efficient physics
simulation codes, especially for modeling strongly correlated electron
systems.
Solution method: We present C++ open source computational libraries that
provide a convenient set of components for developing parallelized physics
simulation codes. The library features a short development cycle and up-to-
date user documentation.

##  1 Introduction

The open source ALPS (Algorithms and Libraries for Physics Simulations)
project [1, 2, 3, 4] provides a collection of computational physics libraries
and applications, especially suited for the simulation of lattice models and
strongly correlated systems. The libraries are lightweight, thoroughly tested
and documented. The libraries are targeted at code developers who implement
generic algorithms and utilities for rapid development of efficient
computational physics applications.

Computer codes based on the ALPS (core) libraries have provided physical
insights into many subfields of condensed matter. Highlights include
nonequilibrium dynamics [5, 6], continuous-time quantum Monte Carlo [7, 8, 9,
10], LDA+DMFT materials simulations [11], simulations of quantum [12] and
classical [13] spin systems, correlated boson [14] and fermion [15] models, as
well as cuprate superconductivity [16].

In this paper we present an updated and refactored version of the ALPS core
libraries. This version introduces a new module to the library: ALPS ALEA,
which is designed as a replacement of the previous accumulator framework. It
provides significant runtime speed-ups for all of its accumulator types,
improved compile times, as well as support for complex random variables, non-
linear propogation of uncertainty and statistical hypothesis testing as a way
to compare stochastic results (Section 2).

Improvements and bug fixes were made to the existing modules. The Green’s
function library was rewritten to be faster and more versatile (Section 3).
Section 4 briefly discusses the other components of the library with the
corresponding improvements.

An effort was made to make the installation procedure more robust and
straight-forward. In Section 5, we list the software packages prerequisite for
compiling and installing the library, as well as outline the installation
procedure. In Section 6 we specify the ALPS citation policy and license, and
in Sections 7 and 8 we summarize the paper and acknowledge our contributors.
Additionally, A provides installation examples.

##  2 ALEA library

The ALEA library (namespace alps::alea) is the next-generation statistics
library for Markov chain Monte Carlo calculations in the ALPS core libraries.
It is intended to supersede the old accumulators library (namespace
alps::accumulators), which is still supported for backwards compatibility.
Compared to the old library, ALEA features:

  * 1.

improved runtime scaling for full binning and logarithmic binning
accumulators, which were both reduced from ${\mathcal{O}}(N\log N)$ to
${\mathcal{O}}(N)$ in terms of the number of samples $N$.[17]

  * 2.

improved runtime performance for vector-valued accumulators due to the use of
the Eigen library and the avoidance of the creation of temporary objects.

  * 3.

improved memory footprint: (a) the use of finalize() semantics avoids
duplicating the accumulator’s data when evaluating the mean and variance; (b)
computed objects allow creating and adding a new measurement in one step
without the need to create a temporary.

  * 4.

improved full binning procedure allows the retention of more fine-grained data
while keeping linear runtime scaling.[17]

  * 5.

improved compilation times and error messages: (a) the interface was greatly
simplified, the metaprogramming logic was removed, and the numerics are now in
the compiled library instead of the headers; (b) the compile-time dependence
on MPI and HDF5 was removed and superseded by clearly specified interfaces.

  * 6.

native support for complex random variables, with user choice of computing
either circular error bars or an error ellipse.

  * 7.

support for computing the full covariance matrix.

  * 8.

support for statistical hypothesis testing [18] as a robust way to compare
Monte Carlo results.

Figure 1: CPU time spent per data point when accumulating $N$ data points of a
random vector with $k=100$ components for the legacy accumulators framework
(red crosses) and the ALEA library (black plusses).

To illustrate the improvements, we compare the CPU time needed per data point
for accumulating the autocorrelation time for a random vector with $k=100$
components between the old accumulators framework and ALEA (Figure 1). One can
see that there is both an offset to the scaling, which comes from avoiding
temporaries, as well as a better overall ${\mathcal{O}}(N)$ vs.
${\mathcal{O}}(N\log N)$ scaling.

### Accumulators

In Monte Carlo simulations, we are interested in obtaining statistical
properties of some random vector $X$ with $k$ components by estimating it from
a sequence of $N$ samples or “measurements” $x=(x_{1},x_{2}\ldots,x_{N})$.
Typically, $N$ is so large that retaining the individual sample points is
impossible or at least impractical. ALEA accumulators provide a workaround by
retaining only certain aggregates of the measurements by performing reductions
over the measurements on the fly, one measurement at a time.

ALEA defines five accumulators, which differ in the stored statistical
estimates and associated runtime and memory cost (cf. Table 1). They represent
commonly used cases ranging from fast and memory efficient to slow but
precise. In all accumulators, size() returns $k$, while count() returns $N$.

ALEA name | old name | time | memory | mean | var | cov | tau | batch
---|---|---|---|---|---|---|---|---
|  |  |  | $\langle x\rangle$ | $\sigma_{x}^{2}$ | $\Sigma_{x}$ | $\tau_{\mathrm{int},x}$ | $\bar{x}_{i}$
mean_acc | MeanAccumulator | $Nk$ | $k$ | ✓ |  |  |  |
var_acc | NoBinningAccumulator | $Nk$ | $k$ | ✓ | ✓ |  |  |
cov_acc | — | $Nk^{2}$ | $k^{2}$ | ✓ | ✓ | ✓ |  |
autocorr_acc | LogBinningAccumulator | $Nk$ | $k\log N$ | ✓ | ✓ |  | ✓ |
batch_acc | FullBinningAccumulator | $Nk$ | $Bk$ | ✓ | ✓ | ✓ |  | ✓
Table 1: Accumulators in ALEA. Asymptotic complexity ${\mathcal{O}}(\cdots)$
for memory and runtime in terms of sample size $N$, number of components of
random vector $k$, and number of bins $B$. Available statistical estimates are
marked with “✓”, where $\langle x\rangle$ is the sample mean, $\sigma_{x}^{2}$
is the sample variance, $\Sigma_{x}$ is the sample covariance matrix,
$\tau_{\mathrm{int},x}$ is the logarithmic binning estimate for the integrated
autocorrelation time, and $\bar{x}_{i}$ are bin means.

### Complex variances

In the case of complex random variables, the concept of error bars becomes
somewhat ambiguous. ALEA supports two common strategies by admitting an
additional Str template argument to var_acc and cov_acc:

  * 1.

Error with circularity constraint (circular_var): error based on the distance
from the mean (circle) in the complex plane. Provides an upper bound to the
total error and ignores any error structure in the real and imaginary part of
the variable.

  * 2.

Error ellipse (elliptic_var): Treats the real and imaginary part of a complex
variable as two separate random variables, thus creating an error ellipse
around the mean in the complex plane. Retains the full error structure in the
complex plane and allows one to plot separate error bars for real and
imaginary part.

### Accumulators and results

Figure 2: Accumulator and result types: finite state machine for accumulator
lifecycle. The boxes “empty”, “accumulating” and “finalized” represent the
possible states of the accumulator. Calling the method indicated on the edges
takes the accumulator from one state to another.

Each accumulator type has a matching result type (e.g., mean_acc has a
matching mean_result). Accumulators and results are facade types over a common
base type (e.g., mean_data), but differ conceptually and thus have
complementary functionality: accumulators support adding data to them, while
results allow one to perform transformations, reductions, as well as
extracting statistical estimates such as the sample mean.

To obtain a result from an accumulator, the accumulators provide both a
result() and a finalize() method. The result() method creates an intermediate
result, which leaves the accumulator untouched and thus must involve a copy of
the data, while the finalize() method invalidates the accumulator and
repurposes its internal data as simulation result (it replaces the sum with
the mean and the sum of squares with the variance). A finite state machine for
these semantics is shown in Figure 2.

One can not only add a vector of data points to an accumulator, but also a
computed object, which creates and adds the data to the accumulator. The
result and computed semantics taken together allow one to use ALEA with very
large random vectors which may fit in memory only once.

### Reductions and serializations

For parallel codes, all results support reduction (averaging over elements)
through the reduce() method, which takes an instance of the abstract reducer
interface. Depending on the implementation of the reducer, the reduction is
performed over different instances (threads, processes, etc.) using MPI,
OpenMP, ssh etc. without introducing a hard dependency on any of these
libraries.

Similarly, all results support serialization (converting to permanent format)
though the serialize() method and deserialization through the deserialize()
method, which take an instance of the abstract serializer and deserializer
interface, respectively. Depending on the implementation, serialization to
HDF5 as well as support for the HPX parallelization framework [19] is
provided.

### Nonlinear propagation of uncertainty

Transformations of results can be performed using the transform method. The
definition of transform, somewhat simplified, is:

⬇

Result transform(Strategy strategy, Function f, const Result &x);

where strategy is a tag denoting the strategy for error propagation, f is the
function $f$, and x is the stochastic result.

Care has to be taken to correctly propagate the uncertainties if these
functions are non-linear, since neglecting the proper error propagation leads
to biasing both the mean and the variance of the transformed result. Several
strategies exist for performing these transforms, and we summarize the ones
supported by ALEA in Table 2. Similar to the case of the accumulators, there
is a tradeoff between computer time required (in particular where $f$ has to
be evaluated repeatedly) and requirements on the Result type, such as the
presence of binned data.[17]

If $f$ is a linear function, propagation through linear_prop is exact. For
non-linear functions on moderate number of components $k$, we suggest the use
of Jackknife (jackknife_prop), which strikes a good tradeoff between bias
correction and cost. The strategies sampling_prop and bootstrap_prop will be
implemented in a later version.

propagation strategy |  preimage distribution |  mean bias |  variance bias | $N_{f}$ | requires
---|---|---|---|---|---
$\langle x\rangle$ | $\Sigma_{x}$ | $\tau_{x}$ | $\bar{x}_{i}$
no_prop | any | $N^{-1}$ | — | $1$ | ✓ |  |  |
linear_prop | Gauss | $N^{-1}$ | $1$ | $k$ | ✓ | ✓ | ✓ |
sampling_prop | known | $S^{-\frac{1}{2}}$ | $S^{-\frac{1}{2}}$ | $S$ | ✓ | ✓ | ✓ |
jackknife_prop | any | $B^{-2}$ | $B^{-2}$ | $B$ |  |  |  | ✓
bootstrap_prop | any | $S^{-\frac{1}{2}}$ | $S^{-\frac{1}{2}}$ | $S$ |  |  |  | ✓
Table 2: Common strategies for uncertainty propagation through a function
$Y=f(X)$. The preimage distribution is the distribution of $X$. Bias describes
the asymptotic $O(\cdots)$ worst-case behavior of the bias on the mean in
terms of the bin size $B$, the number of measurements $N$ taken, or the number
of samples $S$ used for resampling. $N_{f}$ summarizes the number of function
invocations required. The statistical estimates required by each strategy are
marked with “✓”. $\langle x\rangle$ is the sample mean, $\Sigma_{x}$ is the
sample covariance matrix, $\tau_{x}^{\mathrm{int}}$ is the integrated
autocorrelation time, and $\bar{x}_{i}$ are the bin means.

### Hypothesis testing

A common task arising in statistical post-processing is to compare a
stochastic result with a deterministic one, e.g., when testing the algorithm
against known benchmark results, as well as comparing two stochastic results
with each other, e.g., when trying to determine whether an iterative procedure
has converged. This is usually done by visual inspection or some _ad hoc_
criterion.

Statistical hypothesis testing[18] provides a controlled alternative to these
methods. ALEA provides the method test_mean(x, y), which allows for hypothesis
testing. If x is a statistical result with proper (co-)variances and y is a
benchmark (or vice-versa), test_mean performs Hotelling’s $T^{2}$ test for the
following null hypothesis $H_{0}$:[20]

$H_{0}:\langle x\rangle=y,$ (1)

whereas if both are statistical results, the null hypothesis

$H_{0}:\langle x\rangle=\langle y\rangle$ (2)

is checked. In both cases, an object providing the test score and the
corresponding $p$-value is returned, which, roughly speaking, indicates the
probability of the two results being in agreement. The $p$-value can then be
compared with a suitable threshold $\alpha\in[0,1)$, and one considers the
results equal (accepts the null hypothesis) if $p>\alpha$.

##  3 Green’s function library

The _Green’s functions_ component provides a type-safe interface to manipulate
objects representing bosonic or fermionic many-body Green’s functions, self-
energies, susceptibilities, polarization functions, and similar objects. From
a programmer’s perspective, these objects are multidimensional arrays of
floating-point or complex numbers, defined over a set of meshes and
addressable by a tuple of indices, each belonging to a grid. Currently, real
frequency, Matsubara (imaginary frequency), imaginary time (uniform, power,
[21] Legendre [22] and Chebyshev [23]), momentum space, real space, and
arbitrary index meshes are supported.

These many-body objects often need to be supplemented with analytic tail
information encapsulating the high frequency / short time moments of the
Green’s functions, so that high precision Fourier transforms, density
evaluations, or energy evaluations can be performed. In addition to this
functionality, the Green’s function component supports saving data to and
loading it from binary HDF5 files. The UML diagram illustrating the Green’s
function classes and their interdependencies is given in Fig. 3.

Figure 3: UML class diagram of Green’s function library, indicating the main
classes in the Green’s function library. Each box corresponds to a class, and
arrows correspond to relationships between classes.

The central class of the Green’s function library is gf_base. This class
should be parametrized by the scalar type, the type of storage and the list of
meshes (grids). It provides the following basic arithmetic operations:

  * 1.

Addition/subtraction;

  * 2.

In-place addition/subtraction;

  * 3.

Multiplication by a constant.

The addition/subtraction operations can be performed on the Green’s function
objects that are defined on the same grid and have the same shape of the data
array. There are two possible storage types for this class. The memory storage
version of the class is called greenf, whereas the view-based class is called
greenf_view. Instead of having only pre-defined Green’s function classes of
fixed dimension, we provide a flexible interface that makes it possible to
declare Green’s function types of any dimension. The following example defines
4-dimensional Green’s function on Matsubara, arbitrary index, imaginary time
and Legendre meshes:

⬇

alps::gf::matsubara_positive_mesh x(100, 10);

alps::gf::index_mesh y(10);

alps::gf::itime_mesh z(100, 10);

alps::gf::legendre_mesh w(100, 10);

greenf<double, alps::gf::matsubara_positive_mesh, alps::gf::index_mesh,
alps::gf::itime_mesh, alps::gf::legendre_mesh> g(x, y, z, w);

Using view-storage Green’s function one can easily get access to a part of
another Green’s function object using index slicing. This idea is illustrated
by the following code,

⬇

alps::gf::matsubara_positive_mesh x(100, 10);

alps::gf::index_mesh y(20);

greenf<double, alps::gf::matsubara_positive_mesh, alps::gf::index_mesh>
g(x,y);

// loop over leading index

for(alps::gf::matsubara_positive_mesh::index_type w(0); w<x.extent(); ++w) {

greenf_view<double, alps::gf::index_mesh> g2 = g(w);

// do some operation on the Green’s function view object

// defined on the ’y’ mesh alone

}

In the current implementation we support only a slicing over a set of the
leading indices. More details on the Green’s function library can be found in
the respective tutorial.

##  4 Some of the other key library components

In this section, we provide a brief overview of the purpose and functionality
of some of the key components comprising the ALPS core libraries. For a more
detailed exposition with a working example, see the previous publication[1].
Additional information is available online and as part of the Doxygen code
documentation. It should be noted that the components maintain minimal
interdependence, and using a component in a program does not bring in other
components, unless they are required dependencies, in which case they will be
used automatically.

### Building a parallel Monte Carlo simulation

A generic Monte Carlo simulation can be easily assembled from the classes
provided by the _Monte Carlo Scheduler_ component. The programmer needs to
define only the problem-specific methods that are called at each Monte Carlo
step, such as methods to update the configuration of the Markov chain and to
collect the measured data. The simulation is parallelized implicitly, using
one MPI process per chain.

### Storing, restoring, and checkpointing simulation results

To store the results of a simulation in a cross-platform format for subsequent
analysis, one can use the _Archive_ component. The component provides
convenient interface to saving and loading of common C++ data structures
(primitive types, complex numbers, STL vectors and maps), as well as of
objects of user-defined classes to/from HDF5 [24], which is a universally
supported and machine independent data format.

### Reading command-line arguments and parameter files

Input parameters to a simulation can be passed via a combination of a
parameter file and command line arguments. The _Parameters_ library component
is responsible for parsing the files and the command line, and providing
access to the data in the form of an associative array (akin to C++ `map` or
Python dictionary). The parameter files use the standard “`*.ini`” format, a
plain text format with a line-based syntax containing `key = value` pairs,
optionally divided into sections.

##  5 Prerequisites and Installation

To build the ALPS core libraries, any recent C++ compiler can be used; the
libraries are tested with GCC [25] 4.8.1 and above, Intel [26] C++ 15.0 and
above, and Clang [27] 3.2 and above. The library follows the C++11 standard
[28] to facilitate the portability to a wide range of programming
environments, including HPC clusters with older compilers. The library depends
on the following external packages:

  * 1.

The CMake build system [29] of version 3.1 and above.

  * 2.

The Boost C++ libraries [30] of version 1.56.0 and above. Only the headers of
the Boost library are required.

  * 3.

The HDF5 library [24] version 1.8 and above.

  * 4.

The Eigen library [31] version 3.3.4 and above (can be requested to be
downloaded automatically).

To make use of (optional) parallel capabilities, an MPI implementation
supporting standard 2.1 [32] and above is required. Generating the developer’s
documentation requires Doxygen [33] along with its dependencies.

The installation of the ALPS core libraries follows the standard procedure for
any CMake-based package. The first step is to download the ALPS core libraries
source code; the recommended way is to download the latest ALPS core libraries
release from https://github.com/ALPSCore/ALPSCore/releases. Assuming that all
above-mentioned prerequisite software is installed, the installation consists
of unpacking the release archive and running CMake from a temporary build
directory, as outlined in the shell session example below (the `$` sign
designates a shell prompt):

⬇

$ tar -xzf ALPSCore-2.3.0.ta#()r.gz

$ mkdir build

$ cd build

$ export ALPSCore_DIR=$HOME/software/ALPSCore

$ cmake -DCMAKE_INSTALL_PREFIX=$ALPSCore_DIR \

-DALPS_INSTALL_EIGEN=yes \

../ALPSCore-2.3.0

$ make

$ make t#()est

$ make install

The command at line 1 unpacks the release archive (version 2.3.0 in this
example); at line 4 the destination install directory of the ALPS core
libraries is set (`$HOME/software/ALPSCore` in this example). At line 6 the
downloading and co-installation of Eigen library is requested (see A for more
details).

The ALPS core libraries come with an extensive set of tests; it is strongly
recommended to run the tests (via `make test`) to verify the correctness of
the build, as it is done at line 9 in the example above.

The installation procedure is outlined in more details in A; Also, the file
`common/build/build.jenkins.sh` in the library release source tree contains a
build and installation script that can be further consulted for various build
options.

Binary packages are available for some operating systems. On macOS operating
system, the ALPS core libraries package can be downloaded and installed from
the MacPorts [34] repository, using a command `port install alpscore`. On
GNU/Debian Linux operating system, the ALPS core Debian package is provided by
the MateriApps LIVE! project [35] (see [36] for more details).

##  6 License and citation policy

The GitHub version of ALPS core libraries is licensed under the GNU General
Public License version 2 (GPL v. 2) [37] or later. The older ALPS license
under which previous versions of the code were licensed [4] has been retired.
We kindly request that the present paper be cited, along with any relevant
original physics or algorithmic paper, in any published work utilizing an
application code that uses this library.

##  7 Summary

We have presented an updated and repackaged version of the core ALPS
libraries, a lightweight C++ library, designed to facilitate rapid development
of computational physics applications, and have described its main features.

The new version contains the next-generation statistics library ALEA, a major
overhaul of the Green’s function library as well as updates to other
components.

##  8 Acknowledgments

Work on the ALPS library project is supported by the Simons collaboration on
the many-electron problem. Aspects of the library are supported by NSF
DMR-1606348 and DOE ER 46932.

## References

  * [1] A. Gaenko, A. Antipov, G. Carcassi, T. Chen, X. Chen, Q. Dong, L. Gamper, J. Gukelberger, R. Igarashi, S. Iskakov, M. Könz, J. LeBlanc, R. Levy, P. Ma, J. Paki, H. Shinaoka, S. Todo, M. Troyer, E. Gull, Updated core libraries of the alps project, Comput. Phys. Commun. 213 (2017) 235 – 251\.  doi:10.1016/j.cpc.2016.12.009.
URL http://www.sciencedirect.com/science/article/pii/S0010465516303885

  * [2] F. Alet, P. Dayal, A. Grzesik, A. Honecker, M. Körner, A. Läuchli, S. R. Manmana, I. P. McCulloch, F. Michel, R. M. Noack, G. Schmid, U. Schollwöck, F. Stöckli, S. Todo, S. Trebst, M. Troyer, P. Werner, S. Wessel, The ALPS project: Open source software for strongly correlated systems, Journal of the Physical Society of Japan 74 (Suppl) (2005) 30–35.  arXiv:cond-mat/0410407, doi:10.1143/JPSJS.74S.30.
URL http://dx.doi.org/10.1143/JPSJS.74S.30

  * [3] A. Albuquerque, F. Alet, P. Corboz, P. Dayal, A. Feiguin, S. Fuchs, L. Gamper, E. Gull, S. Gürtler, A. Honecker, R. Igarashi, M. Körner, A. Kozhevnikov, A. Läuchli, S. Manmana, M. Matsumoto, I. McCulloch, F. Michel, R. Noack, G. Pawłowski, L. Pollet, T. Pruschke, U. Schollwöck, S. Todo, S. Trebst, M. Troyer, P. Werner, S. Wessel, The ALPS project release 1.3: Open-source software for strongly correlated systems, Journal of Magnetism and Magnetic Materials 310 (2, Part 2) (2007) 1187 – 1193, proceedings of the 17th International Conference on Magnetism.  doi:10.1016/j.jmmm.2006.10.304.
URL http://www.sciencedirect.com/science/article/pii/S0304885306014983

  * [4] B. Bauer, L. D. Carr, H. G. Evertz, A. Feiguin, J. Freire, S. Fuchs, L. Gamper, J. Gukelberger, E. Gull, S. Guertler, A. Hehn, R. Igarashi, S. V. Isakov, D. Koop, P. N. Ma, P. Mates, H. Matsuo, O. Parcollet, G. Pawłowski, J. D. Picon, L. Pollet, E. Santos, V. W. Scarola, U. Schollwöck, C. Silva, B. Surer, S. Todo, S. Trebst, M. Troyer, M. L. Wall, P. Werner, S. Wessel, The ALPS project release 2.0: open source software for strongly correlated systems, Journal of Statistical Mechanics: Theory and Experiment 2011 (05) (2011) P05001.
URL http://stacks.iop.org/1742-5468/2011/i=05/a=P05001

  * [5] M. Eckstein, M. Kollar, P. Werner, Thermalization after an interaction quench in the Hubbard model, Phys. Rev. Lett. 103 (2009) 056403.  doi:10.1103/PhysRevLett.103.056403.
URL http://link.aps.org/doi/10.1103/PhysRevLett.103.056403

  * [6] P. Werner, T. Oka, A. J. Millis, Diagrammatic Monte Carlo simulation of nonequilibrium systems, Phys. Rev. B 79 (2009) 035320.  doi:10.1103/PhysRevB.79.035320.
URL http://link.aps.org/doi/10.1103/PhysRevB.79.035320

  * [7] P. Werner, A. Comanac, L. de’ Medici, M. Troyer, A. J. Millis, Continuous-time solver for quantum impurity models, Phys. Rev. Lett. 97 (2006) 076405.  doi:10.1103/PhysRevLett.97.076405.
URL http://link.aps.org/doi/10.1103/PhysRevLett.97.076405

  * [8] E. Gull, P. Werner, O. Parcollet, M. Troyer, Continuous-time auxiliary-field Monte Carlo for quantum impurity models, EPL (Europhysics Letters) 82 (5) (2008) 57003.
URL http://stacks.iop.org/0295-5075/82/i=5/a=57003

  * [9] H. Shinaoka, E. Gull, P. Werner, Continuous-time hybridization expansion quantum impurity solver for multi-orbital systems with complex hybridizations, Computer Physics Communications 215 (2017) 128–136.  doi:10.1016/j.cpc.2017.01.003.
URL http://linkinghub.elsevier.com/retrieve/pii/S0010465517300036

  * [10] H. Shinaoka, Y. Nomura, E. Gull, Efficient implementation of the continuous-time interaction-expansion quantum monte carlo method.  arXiv:1807.05238.
  * [11] J. Ferber, K. Foyevtsova, R. Valentí, H. O. Jeschke, LDA $+$ DMFT study of the effects of correlation in LiFeAs, Phys. Rev. B 85 (2012) 094505\.  doi:10.1103/PhysRevB.85.094505.
URL http://link.aps.org/doi/10.1103/PhysRevB.85.094505

  * [12] E. Cremades, S. Gómez-Coca, D. Aravena, S. Alvarez, E. Ruiz, Theoretical study of exchange coupling in 3d-Gd complexes: Large magnetocaloric effect systems, Journal of the American Chemical Society 134 (25) (2012) 10532–10542, pMID: 22640027\.  doi:10.1021/ja302851n.
URL http://dx.doi.org/10.1021/ja302851n

  * [13] D. Bergman, J. Alicea, E. Gull, S. Trebst, L. Balents, Order-by-disorder and spiral spin-liquid in frustrated diamond-lattice antiferromagnets, Nat Phys 3 (7) (2007) 487–491.
URL http://dx.doi.org/10.1038/nphys622

  * [14] L. Pollet, J. D. Picon, H. P. Büchler, M. Troyer, Supersolid phase with cold polar molecules on a triangular lattice, Phys. Rev. Lett. 104 (2010) 125302.  doi:10.1103/PhysRevLett.104.125302.
URL http://link.aps.org/doi/10.1103/PhysRevLett.104.125302

  * [15] J. P. F. LeBlanc, A. E. Antipov, F. Becca, I. W. Bulik, G. K.-L. Chan, C.-M. Chung, Y. Deng, M. Ferrero, T. M. Henderson, C. A. Jiménez-Hoyos, E. Kozik, X.-W. Liu, A. J. Millis, N. V. Prokof’ev, M. Qin, G. E. Scuseria, H. Shi, B. V. Svistunov, L. F. Tocchio, I. S. Tupitsyn, S. R. White, S. Zhang, B.-X. Zheng, Z. Zhu, E. Gull, Solutions of the two-dimensional Hubbard model: Benchmarks and results from a wide range of numerical algorithms, Phys. Rev. X 5 (2015) 041041.  doi:10.1103/PhysRevX.5.041041.
URL http://link.aps.org/doi/10.1103/PhysRevX.5.041041

  * [16] E. Gull, O. Parcollet, A. J. Millis, Superconductivity and the pseudogap in the two-dimensional Hubbard model, Phys. Rev. Lett. 110 (2013) 216405.  doi:10.1103/PhysRevLett.110.216405.
URL http://link.aps.org/doi/10.1103/PhysRevLett.110.216405

  * [17] M. Wallerberger, E. Gull, Estimating uncertainties in Markov Chain Monte Carlo simulations, submitted to Am. J. Phys.
  * [18] M. Wallerberger, E. Gull, Hypothesis testing of scientific monte carlo calculations, Phys. Rev. E 96 (2017) 053303.  doi:10.1103/PhysRevE.96.053303.
URL https://link.aps.org/doi/10.1103/PhysRevE.96.053303

  * [19] H. Kaiser, T. Heller, B. Adelstein-Lelbach, A. Serio, D. Fey, HPX: A Task Based Programming Model in a Global Address Space, in: Proceedings of the 8th International Conference on Partitioned Global Address Space Programming Models, PGAS ’14, ACM, New York, NY, USA, 2014, pp. 6:1–6:11.  doi:10.1145/2676870.2676883.
URL http://doi.acm.org/10.1145/2676870.2676883

  * [20] H. Hotelling, The generalization of Student’s ratio 2 (3) (1931) 360–378.  doi:10.1214/aoms/1177732979.
URL http://dx.doi.org/10.1214/aoms/1177732979

  * [21] W. Ku, Electronic excitations in metals and semiconductors: Ab initio studies of realistic many-particle systems., Ph.D. thesis, University of Tennessee (2000).
URL http://trace.tennessee.edu/utk_graddiss/2030

  * [22] L. Boehnke, H. Hafermann, M. Ferrero, F. Lechermann, O. Parcollet, Orthogonal polynomial representation of imaginary-time green’s functions, Phys. Rev. B 84 (2011) 075145.  doi:10.1103/PhysRevB.84.075145.
URL https://link.aps.org/doi/10.1103/PhysRevB.84.075145

  * [23] E. Gull, S. Iskakov, I. Krivenko, A. A. Rusakov, D. Zgid, Chebyshev polynomial representation of imaginary-time response functions, Phys. Rev. B 98 (2018) 075127.  doi:10.1103/PhysRevB.98.075127.
URL https://link.aps.org/doi/10.1103/PhysRevB.98.075127

  * [24] The HDF Group.  Hierarchical Data Format, version 5 [online] (1997–2016).  http://www.hdfgroup.org/HDF5/.
  * [25] R. M. Stallman, the GCC Developer Community, Using the GNU Compiler Collection (1988–2016).
URL https://gcc.gnu.org/onlinedocs/gcc-6.1.0/gcc.pdf

  * [26] Intel Corporation, User and Reference Guide for the Intel C++ Compiler 15.0 (2015).
URL https://software.intel.com/en-us/compiler_15.0_ug_c

  * [27] clang: a C language family frontend for LLVM [online] (2016).  http://clang.llvm.org/.
  * [28] ISO, ISO/IEC 14882:2011 Information technology — Programming languages — C++, International Organization for Standardization, Geneva, Switzerland, 2012.
URL
http://www.iso.org/iso/iso_catalogue/catalogue_tc/catalogue_detail.htm?csnumber=50372

  * [29] K. Martin, B. Hoffman, Mastering CMake, Kitware, Inc, 2015.
  * [30] Boost software libraries [online] (2016).  http://www.boost.org/.
  * [31] G. Guennebaud, B. Jacob, et al.  Eigen v3 [online] (2010).  http://eigen.tuxfamily.org.
  * [32] Message Passing Interface Forum, MPI: A Message-Passing Interface Standard, Version 2.1, High Performance Computing Center Stuttgart (HLRS), 2008, https://www.mpi-forum.org/docs/mpi-2.1/mpi21-report.pdf.
  * [33] D. van Heesch.  Doxygen [online] (1997–2015).  http://www.doxygen.org/.
  * [34] The MacPorts Project.  Macports [online] (2007–2014).  http://www.macports.org/.
  * [35] MateriApps LIVE! [online] (2013–2018).  https://cmsi.github.io/MateriAppsLive/.
  * [36] ALPS Core Debian Package [online] (2018).  https://github.com/cmsi/ma-alpscore.
  * [37] GNU general public license.
URL http://www.gnu.org/licenses/gpl.html

##  Appendix A Detailed installation procedure

In the following discussion we assume that all prerequisite software (section
5) is installed, and the ALPS core libraries release (here, release 2.3.0) is
downloaded into the current directory as `ALPSCore-2.3.0.tar.gz`. Also, we
assume that the libraries are to be installed in `ALPSCore` subdirectory of
the current user’s home directory. The commands are given assuming `bash` as a
user shell.

The first step is to unpack the release archive and set the desired install
directory:

⬇

$ tar -xzf ALPSCore-2.3.0.t#()ar.gz

$ export ALPSCore_DIR=$HOME/software/ALPSCore

The next step is to perform the build of the library (note that the build
should not be performed in the source directory):

⬇

$ mkdir build

$ cd build

$ cmake ../ALPSCore-2.3.0 -DCMAKE_INSTALL_PREFIX=$ALPSCore_DIR \

-DALPS_INSTALL_EIGEN=yes

The `cmake` command at lines 5 and 6 accepts additional arguments in the
format -Dvariable=value. A number of relevant CMake variables is listed in
Table 3. The installation process is also affected by environment variables,
some of which are listed in Table 4; the CMake variables take precedence over
the environment variables. The build and installation script
`common/build/build.jenkins.sh` in the ALPS core libraries release source tree
provides an example of using some of the build options.

It should be noted that line 6 of the example above will cause the ALPSCore
installation process to download a copy of Eigen library[31] and co-install it
with ALPSCore. If the Eigen library is already installed on your system, it
may be preferable to use the installed version. In this case, instead of
requesting the installation of Eigen, the location of the Eigen library should
be specified:

⬇

$ cmake ../ALPSCore-2.3.0 -DCMAKE_INSTALL_PREFIX=$ALPSCore_DIR \

-DEIGEN3_INCLUDE_DIR=/usr/l#()ocal/Eigen3

In this example, the location of the Eigen library is assumed to be
/usr/local/Eigen3; the actual location depends on your local Eigen
installation. The directory specified by the -DEIGEN3_INCLUDE_DIR option must
contain the Eigen subdirectory.

Variable | Default value |  Comment
---|---|---
CMAKE_CXX_COMPILER | (system default) |  Path to C++ compiler executable.*
ALPS_CXX_STD | c++11 |  The C++ standard to compile ALPSCore with.
CMAKE_INSTALL_PREFIX | /usr/local |  library target install directory.
CMAKE_BUILD_TYPE | RelWithDebInfo |  Specifies build type.
BOOST_ROOT |  |  Boost install directory. Set if CMake fails to find Boost.
Boost_NO_SYSTEM_PATHS | false |  Set to true to disable search in default system directories, if the wrong version of Boost is found.
Boost_NO_BOOST_CMAKE | false |  Set to true to disable search for Boost CMake file, if the wrong version of Boost is found.
Documentation | ON |  Build developer’s documentation.
ENABLE_MPI | ON |  Enable MPI build (set to OFF to disable).
Testing | ON |  Build unit tests (recommended).
ALPS_BUILD_TYPE | dynamic |  Can be dynamic or static: build libraries as dynamic (“shared”) or static libraries, respectively.*
EIGEN3_INCLUDE_DIR |  |  Location of Eigen headers* (containing Eigen subdirectory)
Eigen3_DIR |  |  Location of CMake-based Eigen installation* (containing Eigen3Config.cmake)
ALPS_INSTALL_EIGEN | NO |  set to yes to request downloading and co-installation of Eigen*
ALPS_EIGEN_UNPACK_DIR | (autodetected) |  directory where to unpack Eigen for co-installation
ALPS_EIGEN_TGZ_FILE | (autodetected) |  location of Eigen archive to unpack
*Note: For the change of this variable to take effect, remove your build directory and redo the build.
Table 3: CMake arguments relevant to building of ALPS core libraries. Variable |  Comment
---|---
CXX |  Path to C++ compiler executable.*
BOOST_ROOT |  Boost install directory. Set if CMake fails to find Boost.
HDF5_ROOT |  HDF5 install directory. Set if CMake fails to find HDF5.
Eigen3_DIR |  Eigen install directory. Set if you use CMake-based installation of Eigen.
*Note: For the change of this variable to take effect, remove your build directory and redo the build.
Table 4: Environment variables arguments relevant to building of ALPS core
libraries.


Control invariant set is critical for guaranteeing safe control and  the problem of computing control invariant set for linear discrete-time system is revisited in this paper by using a data-driven approach. Specifically, sample points on convergent trajectories of linear MPC are recorded, of which the convex hull formulates a control invariant set for the linear system. To approximate the convex hull of multiple sample points, a convex piecewise linear (PWL) fitting framework has been proposed, which yields a polyhedral approximation with predefined complexity. A descent algorithm for the convex PWL fitting problem is also developed, which is guaranteed to converge to a local optimum. The proposed strategy is flexible in computing the control invariant set in high dimension with a predefined complexity. Simulation results show that the proposed data-driven approximation can compute the approximated control invariant set with high accuracy and relatively low computational cost.

§ INTRODUCTION

In scenarios for guaranteeing closed-loop stability of model predictive control (MPC) and safety-critical situations, the calculation of control invariant set is critical <cit.>. A natural thought for guaranteeing safety or closed-loop stability is to drive the system states into the control invariant set, and the states are  guaranteed to be within the control invariant set since then. In recent years, this topic has gained considerable interests, especially in the area of safe learning for reinforcement learning <cit.>.

Traditional method for computing the maximal control invariant set is based on geometric approach, in which the control invariant set is computed iteratively through the one-step backward operator until a termination condition is met <cit.>.  The initial set can be either  the state constraint set or a control invariant set, corresponding to a non-increasing or non-decreasing sequence of sets, and the maximal control invariant set is approximated from outside or inside, respectively, <cit.>.  In general, this kind of method is often prohibitive for high dimensional or complex systems, and it is difficult to guarantee the termination of the set-inclusion procedure <cit.>.

An alternative procedure is to construct a more conservative but more computationally affordable sets, i.e., the control invariant set is assumed to be a 0-symmetric ellipsoid or polytope, and calculating the control invariant set then reduces to solve an optimization problem to maximize the volume of the ellipsoid or polytope. It has been proved in <cit.> that for linear systems with polyhedral constraints, the maximal control invariant set is also polyhedral, hence polyhedral approximation of the maximal control invariant set is more appropriate. Methods for calculating the polyhedral approximations fall into two categories. The first is by generating a set of non-decreasing control invariant set and introducing more flexible terminate conditions <cit.>. The second kind  is the data-based strategy, in which the safe set (actually control invariant) is constructed based on the closed-loop trajectory of the linear system. It was shown in <cit.> that the convex hull of states on a trajectory that terminates at the origin is control invariant. This idea is used in several work to construct safe sets or control invariant set, like  <cit.>, etc. As a matter of fact, when the number of points on the trajectories is large, it is difficult to calculate the convex hull. In <cit.>, a semi-definite programming problem was formulated to solve an approximated ellipsoid with maximal volume, which is computational expensive when the number of sample points is large.

In this paper, based on the idea of <cit.>, we propose a data-driven approach to approximate the maximal control invariant set. For linear systems with 0-symmetric constraints, we prove that 0-symmetric control invariant set exist. Then the convex hull of 0-symmetric convergent sample points of linear MPC are used to formulate  a 0-symmetric control invariant set.
To approximate the convex hull of a large number of points in high dimension, we propose a convex piecewise linear (PWL) fitting framework, i.e., fitting the points with a convex PWL function while satisfying linear constraints.  An efficient procedure for removing redundant constraints is introduced to reduce the complexity of the convex PWL fitting problem and a descent algorithm  is put forward, which can converge to a local optimum of the optimization problem.

The rest of the paper is organized as follows. Section II describes some preliminaries, then sampling on the trajectories of linear MPC problem as well as the formulated control invariant set is introduced in Section III. Section IV describes the main work, i.e., polyhedral approximation of the control invariant set, in which the 0-symmetry of the control invariant set under symmetric constraints is proved, and the convex PWL fitting problem is proposed and solved. Simulation results are given in Section V and the paper ends with conclusions and future work in Section VI.

§ PRELIMINARIES

Some basic definitions and assumptions are needed before proceeding to the computation of the control invariant set for linear systems. In this manuscript, we consider the constrained linear discrete-time system

\begin{equation}\label{eq:LTI}
\bm x_{k+1}=A\bm x_k+B\bm u_k, \bm x \in \mathcal{X}, \bm u \in \mathcal{U},
\end{equation}

in which $\bm x_k \in \mathbb{R}^{n_x}, \bm u_k \in \mathbb{R}^{n_u}$,  the set $\mathcal{X}$ and $\mathcal{U}$ are polyhedral and 0-symmetric. The definition
of 0-symmetric convex set is given in Definition  <ref>.

A  convex set $\mathcal{S}$ is 0-symmetric if

\[
\bm x \in \mathcal{S} \Rightarrow -\bm x \in \mathcal{S}.
\]

The control invariant set for the constrained linear system (<ref>) is defined as follows <cit.>.

The non-empty set $\Omega \subset \mathbb{R}^n$ is a control invariant set for the system (<ref>) if for all $\bm x_k \in \Omega$, there is $\bm u_k \in \mathcal{U}$, such that $\bm x_{k+1} \in \Omega$. The set $\Omega_{\max}$ is said to be maximal control invariant set for the system (<ref>) if $\Omega \subset \Omega_{\max}, \forall \Omega$.

§ SAMPLING ON THE TRAJECTORIES OF LINEAR MPC TO FORMULATE A CONTROL INVARIANT SET

According to <cit.>, the convex hull of points on convergent trajectories constitutes a control invariant set, and here we use linear MPC to generate such convergent trajectories.

§.§ Linear MPC problem and the property of solutions

We first review the formulation of linear MPC problem and discuss the property of the solution.
The MPC problem for linear discrete-time time-invariant systems can be described as follows,

\begin{align}
\min_U~&\bm x_{N_p}^TP\bm x_{N_p}+\sum_{k=0}^{N_p-1}(\bm x_k^TQ\bm x_k+\bm u_k^TR\bm u_k) \label{eq:obj_linearMPC}\\
s.t.~&\bm x_{k+1}=A\bm x_k+B\bm u_k, k=0, 1,\ldots,N_p-1, \label{eq:linear_dynamic}\\
{}&\bm x_k \in \mathcal{X}, k=1,\ldots,N_p, \label{eq:x_constraint}\\
{}&\bm u_k \in \mathcal{U}, k=0, 1,\ldots, N_p-1, \label{eq:u_constraint}\\
{}&\bm x_{N_p} \in \mathcal{X}_f,\label{eq:xn_constraint}
\end{align}

in which $\bm x_0\in \mathbb{R}^{n_x}$ is the initial state, and $N_p$ is the prediction horizon. The optimized variable is $U=[\bm u_0^T, \ldots, \bm u_{N_p-1}^T]^T\in \mathbb{R}^{n_u\cdot N_p}$ and $P, ~Q \in \mathbb{R}^{n_x \times n_x}, R \in \mathbb{R}^{n_u \times n_u}$ denote the weight matrices. The terminal set $\mathcal{X}_f$ is control invariant, i.e., once the state enters $\mathcal{X}_f$, there are feasible inputs such that the subsequent state remains in $\mathcal{X}_f$.

It was shown in <cit.> that the optimal solution $U^*$ is basically a continuous piecewise linear (PWL) function of $\bm x$, i.e., in subregions, $U^*$ admits an affine expression of $\bm x$, then the closed-loop system becomes a continuous PWL system, i.e.,

\begin{equation}\label{eq:closed_loop}
\bm x_{k+1}=A\bm x_k+Bf_{\rm CPWL}(\bm x_k),
\end{equation}

in which $f_{\rm CPWL}$ is a continuous PWL function.

It was shown in <cit.> that if the terminal cost $\bm x_{N_p}^TP\bm x_{N_p}$ and the terminal set $\mathcal{X}_f$ are chosen suitably, the resulting closed-loop system (<ref>) is exponentially stable within a region of attraction. It has also been proved in <cit.> that for stabilizing $(A, B)$, by increasing the prediction horizon $N_p$, one can obtain an optimal solution $U^*$ that can stabilize the linear system (<ref>), by solving the linear MPC problem (<ref>) with the constraint (<ref>) removed. Hence in this paper, we choose a suitable $N_p$ and omit the constraint (<ref>) to generate convergent trajectories, which is more computationally affordable as calculating the terminal invariant set requires additional work.

§.§ Sample points collection

In this subsection, initial sample points are uniformly generated in $\mathcal{X}$, according to Section <ref>, we can choose a long enough prediction horizon $N_p$ such that if the initial point is feasible, then the continuous PWL system (<ref>) is stable and we record all the trajectories that converge to the equilibrium point $\bm 0$.
Assume all the convergent points are,

\[
\mathcal{X}_s=\{\bm x_1,  \ldots, \bm x_{N_s}\},
\]

similar to <cit.>,  we show that the convex hull of $\mathcal{X}_s$ is control invariant for the original linear system (<ref>).

Suppose for the closed-loop system of linear MPC (<ref>), every trajectory starting from $\bm x_i \in \mathcal{X}_s, i=1, \ldots, N_s$ converges to $\bm 0$, then the convex hull of $\bm x_1, \ldots, \bm x_{N_s}$, denoted as $\mathrm{conv}(\mathcal{X}_s)$ is control invariant for the linear system (<ref>). Besides, if the vertices of the  polyhedral maximal control invariant set are all included in $\mathcal{X}_s$, the constructed convex hull $\mathrm{conv}(\mathcal{X}_s)$ is the maximal control invariant set for the linear system (<ref>).

For every $\bm x_i \in \mathcal{X}_s$, there is a control input $\bm u_i$, which is optimal for the MPC problem (<ref>), such that

\[
A\bm x_i+B\bm u_i \in \mathrm{conv}(\mathcal{X}_s).
\]

Then for any $\bm x \in \mathrm{conv}(\mathcal{X}_s)$, i.e., there is a $\lambda$ such that

\begin{equation}\label{eq:cvx_comb}
\bm x=\sum_{i=1}^{N_s}\lambda_i \bm x_i,
\end{equation}

we can set $\bm u(\bm x)=\sum_{i=1}^{N_s}\lambda_i \bm u_i$, then we have

\[
A\bm x+B\bm u \in \mathrm{conv}(\mathcal{X}_s).
\]

Thus $\mathrm{conv}(\mathcal{X}_s)$ is control invariant.

If all the vertices of the polyhedral maximal control invariant set, say $\Omega_{\max}$, are included in $\mathcal{X}_s$, then from the definition of convex hull, we have

\[
\Omega_{\max} \subset \mathrm{conv}(\mathcal{X}_s).
\]

According to the definition of maximal control invariant set, we also have $\mathrm{conv}(\mathcal{X}_s) \subset \Omega_{
\max}$, which directly yields the conclusion that $\Omega_{\max}=\mathrm{conv}(\mathcal{X}_s)$.

Then we demonstrate that a 0-symmetric controlled invariant set can always be obtained by 0-symmetric sample points.

§.§ 0-symmetric control invariant set

We can construct a  0-symmetric control invariant set $\Omega$, i.e., $\forall \bm x \in \Omega$, we have

\[
-\bm x \in \Omega,
\]

if the constraints set $\mathcal{X}$ and $\mathcal{U}$ are 0-symmetric.

According to the definition in <cit.>, if we choose two matrices $\Theta =-I_{n_x\times n_x}$ and $\Gamma =- I_{n_u \times n_u}$, in which $I_{n_x \times n_x}$ denotes the $n_x \times n_x$ identity matrix, then $(\Theta, \Gamma)$ is a symmetry of the dynamics (<ref>), constraints (<ref>)-(<ref>) and the cost (<ref>). Assume the set $\bar{\Omega}\subset \mathcal{X}$ is control invariant, i.e., for any $\bm x \in \bar{\Omega}$, there is some $\bm u \in \mathcal{U}$ such that $\bm x_{\rm next}=A\bm x+B\bm u \in \bar{\Omega}$.

According to <cit.>, as $\bm x$ is feasible, $-\bm x$ is also feasible. Denote $\tilde{\Omega}=\{-\bm x|\bm x \in \bar{\Omega}\}$, then we have

\[
-A\bm x-B\bm u=-\bm x_{\rm next} \in \tilde{\Omega}, \forall \bm x \in \bar{\Omega}.
\]

Let $\Omega = \bar{\Omega}\cup \tilde{\Omega}$, apparently $\Omega$ is 0-symmetric and for any $\bm x\in \Omega$ we have

\[
A\bm x+B\bm u \in \Omega,
\]

indicating that $\Omega$ is control invariant, and the conclusion follows.

According to Lemma <ref> and <ref>, if we have obtained the sample point set $\mathcal{X}_s$, we can augment the sample point set to obtain 0-symmetric sample point set, i.e.,

\begin{equation}\label{eq:sample_aug}
\tilde{\mathcal{X}}_s=\mathcal{X}_s \cup (-\mathcal{X}_s),
\end{equation}

To avoid the abuse of notation, hereafter we use $\mathcal{X}_s$ to represent $\tilde{\mathcal{X}}_s$. For any $\bm x\in \mathrm{conv}(\mathcal{X}_s)$, there is a $\lambda$ such that (<ref>) holds. Then we have

\[
-\bm x=\sum\limits_{i=1}^{N_s}\lambda_i (-\bm x_i),
\]

and $-\bm x  \in \mathcal{X}_s$, i.e., a 0-symmetric control invariant set can be obtained through the convex hull of ${\mathcal{X}}_s$.

§ POLYHEDRAL APPROXIMATION OF THE CONTROL INVARIANT SET

Different from the existing research, which constructs a maximum-volume ellipsoid <cit.>, here by investigating on the  0-symmetry of the control invariant set, we propose a method based on convex PWL fitting, i.e., fitting the convex hull of sample points, i.e., $\mathrm{conv}(\mathcal{X}_s)$, with a convex PWL function.

§.§ Convex PWL fitting problem formulation

As $\Omega$ is 0-symmetric, we only need to construct a polyhedral approximation of a part of $\Omega$. Specifically, we partition $\Omega$ into 3 sets, i.e.,

\begin{equation}\label{eq:partition}
\Omega=\Omega_0 \cup \Omega_{\rm N} \cup \Omega_{\rm P},
\end{equation}

in which the regions are defined as $\Omega_0=\Omega \cap \{\bm x|x_{n_x}=0\}$, and

\[
\begin{split}
\Omega_{\rm N} = \Omega \cap \{\bm x|x_{n_x}<0\},
\end{split}
\]

\[
\begin{split}
\Omega_{\rm P} = \Omega \cap \{\bm x|x_{n_x}>0 \},
\end{split}
\]

respectively. Here $x_j$ denotes the $j$-th component of $\bm x$, $j=1, \ldots, n_x$.

According to the partition (<ref>), all the sample points obtained in Section <ref> are partitioned into 3 parts, i.e., we have index sets $\mathcal{I}_0, \mathcal{I}_{\rm N}$ and $\mathcal{I}_{\rm P}$, which are defined as,

\[
\mathcal{I}_0=\{i|\bm x_i \in \Omega_0\}, \mathcal{I}_{\rm N}=\{i|\bm x_i \in \Omega_{\rm N}\}, \mathcal{I}_{\rm P}=\{i|\bm x_i \in \Omega_{\rm P}\},
\]

and we have

\[
\{\bm x_i|i \in \mathcal{I}_0\} \cup \{\bm x_i|i \in \mathcal{I}_{\rm N}\} \cup \{\bm x_i|i \in \mathcal{I}_{\rm P}\}=\mathcal{X}_s.
\]

As is discussed before, the control invariant set $\Omega$ is 0-symmetric, a convex PWL function is constructed to formulate a lower bound for the points in $\Omega_0 \cup \Omega_{\rm N}$, the boundaries for points in $\Omega_0\cup \Omega_{\rm P}$ is then obtained by the property of symmetry. The detailed procedure is as follows.

For the sample points $\bm x_i, i \in \mathcal{I}_0 \cup \mathcal{I}_{\rm N}$,
a convex PWL fitting is constructed through designing a convex PWL function to approximate the lower bound of the convex hull $\mathrm{conv}(\bm x_i), i \in \mathcal{I}_0 \cup \mathcal{I}_{\rm N}$, i.e.,

\[\hat{x}_{n_x}=\bigvee\limits_{k \in \{1, \ldots, M\}} \ell_k(x_1, \ldots, x_{n_x-1}),
\]

in which $\ell_k(x_1, \ldots, x_{n_x-1})$ is an affine function, such that

\begin{equation}\label{eq:cvx_pwl}
x_{i,n_x}\geq \hat{x}_{i, n_x}, \forall i \in \mathcal{I}_0 \cup \mathcal{I}_{\rm N}.
\end{equation}

Here $x_{i, j}$ is the $j$-th component for sample point $\bm x_i$, $i=1,\ldots, N_s$, $j = 1, \ldots, n_x$, and $\hat{x}_{i, j}$ can be expressed as follows,

\begin{equation}\label{eq:hat_i}
\hat{x}_{i, n_x}=\bigvee_{k \in \{1, \ldots, M\}}\ell_k(x_{i,1}, \ldots, x_{i, n_x-1}).
\end{equation}

Then we have

\begin{equation}\label{eq:inequality1}
\begin{split}
\ell_k(x_{i,1}, \ldots, x_{i, n_x-1})\leq x_{i,n_x}, \forall k \in \{1, \ldots, M\},
\forall i \in \mathcal{I}_0 \cup \mathcal{I}_{\rm N}.
\end{split}
\end{equation}

As for all $\bm x \in \Omega_0 \cup \Omega_{\rm P}$, the coordinate $x_{n_x}$ is larger than that for points $\bm x \in \Omega_0 \cup \Omega_{\rm N}$, hence (<ref>) is satisfied for all $\bm x_i \in \mathcal{X}_s$. Rearranging (<ref>) and we have the lower approximation to $\Omega$, i.e.,

\[
\mathcal{P}\bm x \leq q,
\]

where $\mathcal{P}$ and $q$ are obtained through (<ref>).
According to the 0-symmetry of the control invariant set $\Omega$, we further have

\begin{equation}\label{eq:linear_inequality}
-\mathcal{P} \bm x \leq q,
\end{equation}

and the control invariant set $\Omega$ can be expressed as

\[
\left[\begin{array}{c}
\mathcal{P}\\
-\mathcal{P}\end{array}\right]\bm x\leq \left[\begin{array}{c}
\]

The polyhedral approximation can also be constructed with respect to the points in $\Omega_0 \cup \Omega_{\rm P}$, i.e., $\hat{x}_{n_x}$ can also be constructed as an upper bound for the points with a  nonnegative $x_{n_x}$, i.e.,

\[
\hat{x}_{n_x}= \bigwedge \limits_{k \in \{1, \ldots, M\}} \ell_k(x_{1}, \ldots, x_{n_x-1}), \forall i \in \mathcal{I}_0 \cup \mathcal{I}_{\rm P},
\]

and the points $\bm x_i$ with $i \in \mathcal{I}_0 \cup \mathcal{I}_{\rm P}$ satisfies,

\[
x_{i,n_x}\leq \bigwedge \limits_{k \in \{1, \ldots, M\}} \ell_k(x_{i,1}, \ldots, x_{i,n_x-1}), \forall i \in \mathcal{I}_0 \cup \mathcal{I}_{\rm P}.
\]

This can also been used to get a linear inequality expression, which is actually (<ref>).

The construction of the convex PWL fit to the sample points can be summarized as the following optimization problem, in which the distance between convex PWL function and the sample $x_{i, n_x}$ is minimized subject to the constraints that the sample points are above the convex PWL function.

\begin{equation}\label{eq:cvx_fit_prob}
\begin{array}{rl}
\min\limits_{\ell_1, \ldots, \ell_M} &J=\sum_{i \in \mathcal{I}_0 \cup \mathcal{I}_{\rm N}}\left(\hat{x}_{i, n_x}-x_{i,n_x}\right)^2 \\
s.t. & (\ref{eq:hat_i})~ \mbox{and}~ (\ref{eq:inequality1}). %\label{eq:pwq}
\end{array}
\end{equation}

It is clear that for the above optimization problem,  the constraints are linear and the cost function is piecewise quadratic, i.e.,

\begin{equation}\label{eq:Jquadratic}
J=\sum\limits_{i \in \mathcal{I}_0 \cup \mathcal{I}_{\rm N}}\left(\ell_{k_i}(x_{i, 1}, \ldots, x_{i, n_x-1})-x_{i, n_x}\right)^2,
\end{equation}

for all $\ell_1, \ldots, \ell_M$ such that

\begin{equation}\label{eq:linear_region}
\begin{array}{r}
\ell_{k_i}(x_{i,1}, \ldots, x_{i, n_x-1})\geq \ell_{k}(x_{i,1}, \ldots, x_{i, n_x-1}), \\
\quad\quad \forall i \in \mathcal{I}_0 \cup \mathcal{I}_{\rm N}, \forall k \in \{1, \ldots, M\}, k \neq k_i.
\end{array}
\end{equation}

§.§ Reducing the complexity of the convex PWL fitting problem

As can be seen from the optimization problem (<ref>), the number of linear constraints is $M\cdot(|\mathcal{I}_{\rm N}|+|\mathcal{I}_0|)$, in which $|\cdot|$ denotes the number of elements in a set. When $|\mathcal{I}_{\rm N}|+|\mathcal{I}_0|$ is large, it is difficult to solve the optimization problem (<ref>). Hence we perform a procedure to reduce the number of constraints and the following lemma provides the criterion for removing constraints.

If for some sample point $\bm x_{i_0}=[x_{i_0, 1}, \ldots, x_{i_0, n_x}]$, there exist sample points $\bm x_{j}, j\in \mathcal{J}_{i_0}$, where $\mathcal{J}_{i_0}$ is an index set with respect to $\bm x_{i_0}$, such that  $\bm x_{i_0}$ can be seen as the convex combination of $\bm x_{j}, j\in \mathcal{J}_{i_0}$, then if

\[
\ell_k(x_{j,1}, \ldots, x_{j, n_x-1})\leq x_{j,n_x}, \forall k \in \{1, \ldots, M\}, \forall j \in \mathcal{J}_{i_0},
\]

we have

\[
\ell_k(x_{i_0,1}, \ldots, x_{i_0, n_x-1})\leq x_{i_0,n_x}, \forall k \in \{1, \ldots, M\}.
\]

As  $\bm x_{i_0}$ is the convex combination of $\bm x_{j}, j\in \mathcal{J}_{i_0}$, there exist scalars $\lambda_j \in [0, 1], j \in \mathcal{J}_{i_0}$ such that

\[
\bm x_{i_0}=\sum\limits_{j \in \mathcal{J}_{i_0}} \lambda_j \bm x_j.
\]

Then $\forall k \in \{1, \ldots, M\}$, as $\ell_k$ is affine, we have

\[
\ell_k(x_{i_0,1}, \ldots, x_{i_0, n_x-1})=\sum\limits_{j \in \mathcal{J}_{i_0}} \lambda_j \ell_k(x_{j,1}, \ldots, x_{j, n_x-1}).
\]

Therefore, the following inequality holds,

\[
\begin{split}
&\ell_k(x_{i_0,1}, \ldots, x_{i_0, n_x-1})-x_{i_0,n_x}\\
&=\sum\limits_{j \in \mathcal{J}_{i_0}} \lambda_j \left(\ell_k(x_{j,1}, \ldots, x_{j, n_x-1})-x_{j,n_x}\right)\leq 0,
\end{split}
\]

and the conclusion follows.

A simple way to eliminate redundant points is to formulate a simplex, which requires $n_x+1$ affinely independent points $P_1, \ldots, P_{n_x}, P_{n_x+1}$ in $n_x$-th dimension, here affinely independent means that $P_2-P_1, \ldots, P_{n_x+1}-P_1$ are linearly independent.
As the equilibrium point $\bm 0$ is in the centre of the control invariant set, it satisfies the constraint (<ref>). And if the distance between a samples points and $\bm 0$ is large, the sample point is likely to be on the boundary of the control invariant set. Hence we select $\bm 0$ and $n_x$ sample points with large distance to form a simplex, and the constraints corresponding to  points within the simplex can be removed.  The removal procedure can be repeated until the number of candidate points to form a simplex is less than $n_x+1$.

§.§ Descent algorithm for solving the convex PWL fitting problem

Assume $\tilde{\bm x}=[x_1, \ldots, x_{n_x-1}, 1]^T$, $\ell_k(x_1, \ldots, x_{n_x-1})=\bm \alpha_k^T \cdot \tilde{\bm x}$, in which $\alpha_k \in \mathbb{R}^{n_x}$, then for each sample point $\bm x_i$, we have

\begin{equation}\label{eq:lk}
\ell_k(x_{i,1}, \ldots, x_{i, n_x-1})=\bm \alpha_k^T \cdot \tilde{\bm x}_i,
\end{equation}

in which $\tilde{\bm x}_i=[x_{i,1}, \ldots, x_{i, n_x-1}, 1]^T$.

For the convex PWL fitting problem (<ref>), we propose to solve one quadratic programming at a time.
Algorithm <ref> shows the process of solving the convex PWL fitting problem, which is then shown descent and can converge to a local optimum in Lemma <ref>.

Algorithm for solving convex PWL fitting problem (<ref>)

[1]

Optimization problem (<ref>), sample points $\bm x_i, i \in \mathcal{I}_0 \cup \mathcal{I}_{\rm N}$, specified number of linear pieces $M$, a small enough threshold $\epsilon$.

Solution  $\ell_1, \ldots, \ell_M$.

Randomly generate $M$ vectors $\bm \alpha_{1,0}, \ldots, \bm \alpha_{M,0}$, let ${flag}=1$.

For each $i \in \mathcal{I}_0 \cup \mathcal{I}_{\rm N}$, find $k_{i,0}$ such that $\bm \alpha_{k_{i_0}}^T \cdot \tilde{\bm x}_i=\vee_{k \in \{1, \ldots, M\}}\bm \alpha_{k}^T \cdot \tilde{\bm x}_i$.

Solve the following quadratic programming problem

\begin{equation}\label{eq:qp}
\begin{array}{rl}
\min\limits_{\bm \alpha_1, \ldots, \bm \alpha_M} &\tilde{J} = \sum\limits_{i \in \mathcal{I}_0 \cup \mathcal{I}_{\rm N}} \left(\bm \alpha_{k_{i_0}}^T \cdot \tilde{\bm x}_i-x_{i,n_x}\right)^2 \\
s.t. &  (\ref{eq:hat_i})~\mbox{and}~(\ref{eq:inequality1}) %\nonumber
\end{array}
%  \vspace{-0.4cm}
\end{equation}

to obtain optimal $\bm \alpha_1, \ldots, \bm \alpha_M$.

$J(\bm \alpha_1, \ldots, \bm \alpha_M)<J(\bm \alpha_{1,0}, \ldots, \bm \alpha_{M,0})$
Using binary search in the line $[\mathcal{L}_0, \mathcal{L}]$, in which $\mathcal{L}_0=[\bm \alpha_{1,0}^T, \ldots, \bm \alpha_{M,0}^T]^T$ and $\mathcal{L}=[\bm \alpha_{1}^T, \ldots, \bm \alpha_{M}^T]^T$, to search for a point $\hat{\mathcal{L}}=[\hat{\bm \alpha}_{1}^T, \ldots, \hat{\bm \alpha}_{M}^T]^T$ such that

\begin{equation}\label{eq:decrease_alg}
J(\hat{\bm \alpha}_1, \ldots, \hat{\bm \alpha}_M)<J(\bm \alpha_{1,0}, \ldots, \bm \alpha_{M,0}).
\end{equation}

Let $\bm \alpha_1=\hat{\bm \alpha}_1, \ldots, \bm \alpha_M=\hat{\bm \alpha}_M$.
$\|(\bm \alpha_1, \ldots, \bm \alpha_M)-(\bm \alpha_{1,0}, \ldots, \bm \alpha_{M,0})\|\leq \epsilon$
$\bm \alpha_{k, 0}=\bm \alpha_k, k=1, \ldots, M$.
$\ell_1, \ldots, \ell_M$ according to (<ref>).

Algorithm <ref> is a descent algorithm and can converge to a local optimum of the optimization problem (<ref>).

If Line 7 in Algorithm <ref> is evaluated, a solution $(\bm \alpha_1, \ldots, \bm \alpha_M)$ with reduced value $J$ is obtained.

Then we demonstrate that if

\[
J(\bm \alpha_1,\ldots, \bm \alpha_M)>J(\bm \alpha_{1,0}, \ldots, \bm \alpha_{M,0}),
\]

we can find some $\hat{\mathcal{L}}=\theta \mathcal{L}_0+(1-\theta) \mathcal{L}$ such that (<ref>) holds. This is easy to verify. For the line $[\mathcal{L}_0, \mathcal{L}]$, there is some $\hat{\mathcal{L}}$ such that

\[
\hat{\bm \alpha}_{k_{i,0}}^T \cdot \tilde{\bm x}_i=\vee_{k \in \{1, \ldots, M\}}\hat{\bm \alpha}_k^T\cdot \tilde{\bm x}_i, \forall i \in \mathcal{I}_0 \cup \mathcal{I}_{\rm N},
\]

which means that for $\mathcal{L}_0$ and $\mathcal{L}$, $J$ admits the same quadratic function, i.e.,

\begin{align}
&J(\bm \alpha_{1,0},\ldots, \bm \alpha_{M,0})=\tilde{J}(\bm \alpha_{1,0},\ldots, \bm \alpha_{M,0}), \\
&J(\hat{\bm \alpha}_1, \ldots, \bm   \hat{\bm \alpha}_M)=\tilde{J}(\bm  \hat{\bm \alpha}_1, \ldots, \hat{\bm \alpha}_M).
\end{align}

According to the  optimization problem (<ref>), we have

\[
\tilde{J}(\bm \alpha_1, \ldots, \bm \alpha_M)<\tilde{J}(\bm \alpha_{1,0}, \ldots, \bm \alpha_{M,0}).
\]

As the function $\tilde{J}$ is quadratic, and $\hat{\mathcal{L}} \in [\mathcal{L}_0, \mathcal{L}]$, we further have

\[
\tilde{J}(\hat{\bm \alpha}_1, \ldots, \hat{\bm \alpha}_M)<\tilde{J}(\bm \alpha_{1,0}, \ldots, \bm \alpha_{M,0}).
\]

Combining with (<ref>), we have the following conclusion

\[
{J}(\hat{\bm \alpha}_1, \ldots, \hat{\bm \alpha}_M)<{J}({\bm \alpha}_{1,0}, \ldots, \bm {\alpha}_{M,0}),
\]

meaning that Algorithm <ref> is descent and can converge to a local optimum.

In our previous work <cit.>,  a quadratic programming problem is also solved at a time, but within the linear subregion expressed by (<ref>). Compared with the approach proposed in this paper, additional $(M-1)\cdot (|\mathcal{I}_0|+|\mathcal{I}_{\rm N}|)$ linear inequalities are needed in <cit.> to describe the specific linear region when evaluating Line <ref> in Algorithm <ref>, resulting in a much larger computational complexity.

In practice, we can set several candidate number of affine functions  $M$, which relates to the complexity of the approximated polyhedron. Then for each candidate $M$, several starting points are generated, resulting in a series of local optima. The final approximated polyhedron is the one chosen from all the local optima for all the candidate $M$ that yields the minimum $J$. It is noted that here the complexity of the polyhedral approximation can be predefined by tuning $M$.

§ SIMULATION

In this section, we give two examples to illustrate the polyhedral approximation of the control invariant set, one is in two dimension and the other is in four dimension.

§.§ An example in two dimension

Considering the two-dimensional example introduced in <cit.>, the system dynamics can be written as

\[
\begin{array}{rcl}
\end{array}\right]x_k+\left[\begin{array}{cc}
\end{array}\right]u_k.%\\
\end{array}
\]

The  system constraints are $-1 \leq u_k \leq 1$ and $-1 \leq x_{k}\leq 1$.

The parameters of the MPC problem used to generate sample points are  $Q=\mathrm{diag}(1,1)$, $R=5\times 10^3 \mathrm{diag}(1,1)$, $P=10$, and $N_p=10$. Uniformly generate 300 starting points in the domain $[-1, 1]^2$ and the sample points are generated on trajectories, which then yield 7044 samples in total, as Fig. <ref> shows. Among all the 7044 sample points, only 33 are left after removing redundant constraints. Solving the convex PWL fitting problem for points with a nonnegative $x_2$ generates 5 affine functions, and the resulting polyhedral approximation contains 12 linear inequalities. The maximal control invariant set, which is computed through the MPT toolbox 3.0 <cit.>, the polyhedral approximation are shown in red and yellow, respectively. The processes of sampling, removing redundant constraints, and convex fitting cost 19.4477s, 0.97s, and 9.0012s, respectively. All the computations in this paper are implemented through MATLAB 2016b (MathWorks, USA) on an Apple M1 Max computer.

It can be seen from Fig. <ref> that the approximated control invariant set approaches the genuine maximal control invariant set closely.

Approximation of the control invariant set for Example 1.

§.§ A 4-dimensional example

Consider the following 4-dimensional inverted pendulum on cart example <cit.>, for which the system dynamic is,

\[
\begin{array}{cccc}
\end{array}
\right], B=\left[\begin{array}{c}
\end{array}\right].
\]

The prediction horizon is taken to be $N=10$ and the value of the matrices in the cost function is $Q = \rm diag\{2, 2, 2, 2\}$, $R=1$, and $P = 1$. The system constraints are set to be $\|x\|_{\infty}\leq [1, 1.5, 0.35, 1]^T$ and $|u|\leq 1$. To get the convergent trajectories, we uniformly generate 300 starting points and record the system evolution, and then 0-symmetric points are generated as in (<ref>) and we totally obtain 93954 sample points. After removing redundant points using Lemma <ref>, only 547 points are left for the convex PWL fitting problem (<ref>), the result of which is a collection of 10 affine functions. Combined with the system constraints, we obtain an approximated control invariant set in 4-dimension with 26 linear inequalities constraints. Fig. <ref> shows the projection of the approximated control invariant set as well as the sample points. It can be seen from Fig. <ref> that the approximated control invariant set closely approximate the convex hull of  all the sample points.

[Projection onto $x1-x2$ plane]

[Projection onto $x1-x4$ plane]

[Projection onto $x2-x4$ plane]

[Projection onto $x3-x4$ plane]

Projections of the approximated control invariant.

For this example, the sampling process and the constraints removing process takes 171.16s and 186.5s, respectively. Then setting the number of affine functions $M$ to be 10, 14,  18, and solving 100 convex PWL fitting problems with different starting points for each candidate $M$, the whole procedure costs 1070.1s, which is affordable for offline calculation.
It is noted that for this example, the MPT toolbox failed to calculate the exact control invariant set after 24 hours.

§ CONCLUSIONS

In this paper, we have developed a data-driven strategy for approximating the control invariant set of linear constrained system. If the state and input constraints are 0-symmetric, the control invariant set of a linear time-invariant system is proved to be 0-symmetric. Sample points on 0-symmetric convergent trajectories of linear MPC are collected, and a polyhedral approximation to the convex hull of all the sample points is constituted by formulating a convex PWL fitting problem. An algorithm for solving the convex PWL programming is provided, which is descent and can converge to the local minimum. Simulation results show the efficacy of the proposed strategy.

In the future, we will work towards constituting an approximated control invariant set that is  a subset of the maximal control invariant set.