MTSOS: writing barrier functions

In order to compile MTSOS, the user must provide a function for calculating the barrier at a given a point (a_i,b_i,u_i) . Although it is possible to write the barrier funciton in MATLAB, doing so will incur a significant performance penalty, so users are encouraged to write their barrier function directly in C.

A very brief introduction to barrier methods

A very brief introduction to barrier methods will be provided here and we will obscure many of the details to focus on implementation. For more information on barrier functions for interior point methods see Convex Optimization S 11.2 and S 11.3 .

As outlined in Minimum-Time Speed Optimization Over a Fixed Path the constraints will have the form

$(a_i,b_i,u_i) in tilde{mathcal{C}}_{theta_i}, quad i = 1,ldots, n$ .

For the purposes of explanation we will consider constraints of the more specific form

f(a_i,b_i,u_i) leq C(s_i), quad i = 1,ldots,n ,

where $f: mbox{bf R}^{2+r} to mbox{bf R}$ is convex.

In a barrier method, we reformulate our constraints as a cost in the objective by adding I_-(f(a_i,b_i,u_i)-C(s_i)) to the objective, where

$I_(u)=bigg{begin{array}{cc}0&u leq 0 infty&u > 0end{array}$

is the indicator function. If the constraint is satisfied, it will not contribute to the objective at all, while if it is unsatisfied, the objective will be infinite. Therefore, to minimize the objective, the constraints must be satisfied. Unfortunately, the indictaor is not differentiable so we approximate it by

$I_(u) approx -kappalog(-u), quad mathrm{dom} = -mbox{bf R}_{++}$ ,

where kappa > 0 is a parameter that sets the accuracy of the approximation. As kappa decreases, this approximation approaches the indicator function. Although we strongly encourage the use of the log barrier function, MTSOS is written to allow the use of an alternative barrier function if desired. If you have multiple constraints, say k_i of them at evaluation point , then we evaluate them as

$sum_{j=0}^{k_i}-kappa logleft(C_j(s_i)-f_j(a_i,b_i,u_i)right) = phi_i$ .

In the barrier function you must calculate phi_i , nablaphi_i , and nabla^2phi_i as described below.

MATLAB

If the user must write the barrier function in MATLAB, change the compiler to point to barreir_matlab.c instead of barrier_front.c. This will cause the algorithm to look to barrier_mat.m for the dynamics. The name of this file can be changed in barrier_matlab.c.

[H G F] = barrier(S, S_prime, S_dprime, a, b, u,indicator,kappa,variables)

Inputs

`S`	a matrix of points where the barrier function must be evaluated
`S_prime`	a matrix of derivatives of with respect to , evaulated at the same points as `S`
`S_dprime`	a matrix of second derivatives of with respect to evaluated at the same points as `S`
`b`	a vector of the values for ; the term is needed to calculate $dot{q_i}^2$ if it is necessary; in particular $dot{q_i}^2=s^{prime 2}_i b_i$
`a`	a vector of the values of for ; the term is needed to calculate $ddot{q_i}$ if it is necessary; in particular $ddot{q}_i = s_i^prime a_i+s_i^{primeprime} b_i$
`u`	an matrix containing the , $u_i in mbox{bf R}^r$ corresponding to the control inputs
`indicator`	an integer that can take on values between 0 and 3 signifying which terms need to be calculated 0 : Hessian and gradient are required 1 : only the Hessian is required 2 : only the gradient is required 3 : only is required
`kappa`	, used in the calculation of
`variables`	a vector provided by the user at the call to `MTSOS` in the optional parameter `variables`

Outputs

`H`	a length vector consisting of the Hessian of each for provided in column first order and concatenated; if any of the points lie outside the domain of for any , then `H` should be set to `inf` NOTE: variable order is
`G`	a length vector consisting of the gradient of each for concatenated;if any of the points lie outside the domain of for any , then `G` should be set to `inf` NOTE: variable order is
`F`	a length vector consisting of the evaluation of for ; if any of the points lie outside the domain of for any , then G should be set to `inf`

C

int so_barrier(const double* S, const double* S_prime, const double* S_dprime,: const double* b, const double* a, const double* u,
int S_length, int U_size, int State_size,
int indicator, double kappa, const double* variables, int variables_length,
double* H_barrier, double* G_barrier, double* F_barrier)

Inputs

Note: these are all constant.

`double* S`	a length array of points where the barrier function must be evaluated
`double* S_prime`	a length array of the derivative of with respect to , evaluated at the same points as `S`
`double* S_dprime`	a length array of the second derivative of with respect to , evaluated at the same points as `S`
`double* b`	a length array of the values of for ; the term is needed to calculate $dot{q}_i^2$ if it is necessary; in particular $dot{q}_i^2 = s_i^{prime 2} b_i$
`double* a`	a length array of the values of for ;the term is needed to calculate $ddot{q}_i$ if it is necessary: $ddot{q}_i = s_i^prime a_i+ s_i^{primeprime} b_i$
`double* u`	a length array containing the , $u_i in mbox{bf R}^r$ corresponding to the control inputs at each epoch concatenated
`int S_length`	, the number of discretization points
`int U_size`	, the number of control inputs
`int State_size`	, the size of the configuration space
`int indicator`	an integer that can take on values between 0 and 3 signifying which terms need to be calculated 0 : Hessian and gradient are required 1 : only the Hessian is required 2 : only the gradient is required 3 : only is required
`double kappa`	, used in the calculation of
`double* variables`	an array provided by the user at the call to `MTSOS` in the optional parameter `variables`
`int variables_length`	the length of the array `variables`

Outputs

Allocated, but unitialized, arrays are provided for the user to populate. In particular, if an array needs to be zero, the values must be set to zero, do not assume any values in the arrays.

`RETURN int`	`1` : all values are within the domain, no errors `0` : points lie outside the domain of `-1` : an error occurred and the algorithm will abort
`double* H_barrier`	a length array consisting of the Hessian for each for provided in column first order and concatenated; if any of the points lie outside the domain of for any , then the value of `H_barrier` is irrelevant and the function should return 0
`double* G_barrier`	a length array consisting of the gradient of each for concatenated; if any of the points lie outside the domain of for any , then the value of `G_barrier` is irrelevant and the function should return 0
`double* F_barrier`	a length array consisting of the evaluation of for ; if any of the points lie outside the domain of for any , then the value of `F_barrier` is irrelevant and the function should return 0

Tips for writing barrier functions

Log barrier gradients and Hessians

Some useful quantities when using the log barrier function, $phi(x)=-sum_{j=1}^{k}kappalog(-g_j(x)),$ are:

$nablaphi(x) = sum_{j=1}^kfrac{-kappa}{g_j(x)}nabla g_j(x),$

$nabla^2phi(x) = sum_{j=1}^kfrac{kappa}{g_j(x)^2}nabla g_j(x) nabla g_j(x)^T + sum_{i=1}^kfrac{-kappa}{g_j(x)}nabla^2 g_j(x),$

where g_j is f_j(a_i,b_i,u_i)-C_j(s_i) .

This suggests an advantageous structuring of the code for the calculation of the gradient and Hessian of phi .

Example structure

- calculate $X_{(i,j)} = f_{(i,j)} - C_{(i,j)}$ ;
- if ( $X_{(i,j)} > 0$ )
  - return 0;
if (indicator )
- calculate $nabla X_{(i,j)}, quad forall; i,j$ ;
- if (indicator )
  - calculate , known as H
- if (indicator or indicator )
  - calculate , known as G
else
- calculate , known as F