This code is an implementation of the subproblem solver that is presented in the article
Marvin Severitt and Paul Manns. Efficient solution of discrete subproblems arising in integer optimal control with total variation regularization. INFORMS Journal on Computing, 35 (4), 2023.
for subproblems that appear in the sequential linear integer programming algorithm (SLIP) for TV-regularized integer optimal control problems, which is presented in the article
Sven Leyffer and Paul Manns. Sequential linear integer programming for integer optimal control with total variation regularization. ESAIM: Control, Optimisation and Calculus of Variations (ESAIM: COCV) 28 (66), 2022.
Currently, the implementation only supports uniform discretizations of a one-dimensional control domain (a,b)
.
trs4slip
: Main librarytrs4slip_test
: Unit teststrs4slip_python_sdist
: Source package for python bindingstrs4slip_python_whl
: Compiled (wheel) package for python bindings
For the python bindings, you will need a python interpreter with the build
, numpy
, and dev
packages being installed. Currently, python>=3.10 is required. You can also change the version requirement in the CMakeLists.txt to 3.8.
$ mkdir -p build
$ cd build
$ cmake .. -DCMAKE_BUILD_TYPE=Release
$ make
$ sudo make install # <-- system-wide install
Test if the algorithm works correctly after building the C++ library (paths relative to project root)
$ cd build
$ make test
(load your anaconda workspace, etc.)
$ mkdir -p build
$ cd build
$ cmake .. -DCMAKE_BUILD_TYPE=Release
$ make trs4slip_python_whl
$ cd ../python/dist
$ ls
Now, there is a file SOME_COMPLICATED_FILENAME.whl
in this folder. The python module can be installed now
$ pip install SOME_COMPLICATED_FILENAME.whl
$ cd examples
$ python run_slip_algorithm.py
The python package provides a function run
that encapsulates
the C++-function trs4slip_astar
, which hints at the fact that a variant
of the A* algorithm is used to solve the subproblems.
The C++ interface is as follows.
void trs4slip_astar(
int32_t * x_next_out,
double const *c,
int32_t const *x,
int32_t const *bangs,
int const delta,
int const N,
int const M,
double * vert_costs_buffer,
int32_t * vert_layer_buffer,
int32_t * vert_value_buffer,
int32_t * vert_prev_buffer,
int32_t * vert_remcap_buffer,
const int boundcon,
double lbound,
double rbound
);
x_next_out
: used to store (and return) the solution computed by the subproblem solver. The user needs to provide an integer array of size N if N denotes the number of discretization intervals for the control. The quantity corresponds tox + d
in [1].c
: used to store the cost coefficients (see (TR-IP) on p. 4 in [1]). The user needs to provide the cost as a double array of size N. The quantity corresponds toc
in [1].x
: used to store the current iterate. The user needs to provide the current iterate as an integer array of size N if N denotes the number of discretization intervals for the control. The integers need to be elements of bangs (below). The quantity corresponds tox
in [1].bangs
: used to store the feasible control realizations. The user needs to provide an integer array of size M if N denotes the number of feasible control realizations. The quantity corresponds to Ξ in [1].delta
: used to store the trust-region radius as an integer. It is the L1-norm-value scaled by N.N
: number of discretization intervals as integer.M
: number of discrete realizations as integer.vert_costs_buffer
: buffer for internal purposes in order to serve frequent reallocation of memory. The user needs to provide a double array of sizeN*M*(Delta_max + 1) + 2
, whereDelta_max
is the maximum trust-region radius delta that is ever used by the caller.vert_layer_buffer
: buffer for internal purposes in order to serve frequent reallocation of memory. The user needs to provide an integer array of sizeN*M*(Delta_max + 1) + 2
, whereDelta_max
is the maximum trust-region radius delta that is ever used by the caller.vert_value_buffer
: buffer for internal purposes in order to serve frequent reallocation of memory. The user needs to provide an integer array of sizeN*M*(Delta_max + 1) + 2
, whereDelta_max
is the maximum trust-region radius delta that is ever used by the caller.vert_prev_buffer
: buffer for internal purposes in order to serve frequent reallocation of memory. The user needs to provide an integer array of sizeN*M*(Delta_max + 1) + 2
, whereDelta_max
is the maximum trust-region radius delta that is ever used by the caller.vert_remcap_buffer
: buffer for internal purposes in order to serve frequent reallocation of memory. The user needs to provide an integer array of sizeN*M*(Delta_max + 1) + 2
, whereDelta_max
is the maximum trust-region radius delta that is ever used by the caller.boundcon
: flag that indicates if Dirichlet boundary conditions are supplied for x (implying that jumps to these boundary values are also taken into account in the objective value).lbound
: left boundary value (only read if boundcon not false).rbound
: right boundary value (only read if boundcon not false).
The python interface is as follows.
def run(np.int32_t[::1] x_next not None,
double[::1] c not None,
np.int32_t[::1] x not None,
np.int32_t[::1] bangs not None,
int delta,
double[::1] vert_costs_buffer not None,
np.int32_t[::1] vert_layer_buffer not None,
np.int32_t[::1] vert_value_buffer not None,
np.int32_t[::1] vert_prev_buffer not None,
np.int32_t[::1] vert_remcap_buffer not None,
int boundcon,
double lbound,
double rbound
):
The variables are the same as in the C++ case but N and M are deduced from the given numpy arrays. In python, the allocation of the working memory arrays may look as follows.
vert_costs_buffer = np.empty(N*M*(Delta_max + 1) + 2)
vert_layer_buffer = np.empty(N*M*(Delta_max + 1) + 2, dtype=np.int32)
vert_value_buffer = np.empty(N*M*(Delta_max + 1) + 2, dtype=np.int32)
vert_prev_buffer = np.empty(N*M*(Delta_max + 1) + 2, dtype=np.int32)
vert_remcap_buffer = np.empty(N*M*(Delta_max + 1) + 2, dtype=np.int32)