External C/C++ Function Support

The previous sections mainly introduced how to explicitly define problems through symbolic expressions. When the optimization problem cannot be defined symbolically (e.g., when constraints involve neural networks or complex robot dynamics models), it can be defined directly through external C/C++ functions.

Custom Functions

OPTIMake supports C/C++ definitions of the following functions through the external interface:

• objective
• end_objective
• inequality
• equality
• start_equality
• end_equality

When defining external functions, you need to provide the function itself and its Jacobian implementation. In certain scenarios (depending on the hessian_approximation configuration in option), you also need to provide the Hessian implementation:

hessian_approximation Configuration	Does external objective require Hessian	Does external constraint require Hessian
'gauss-newton'	Yes	No
'exact'	Yes	Yes
'bfgs'	No	No

After using the external interface to define an objective or constraint as an external function, the code generation feature will automatically generate a C/C++ template for that external function.

Below is the template for an objective's external C/C++ function:

void vehicle_l(
    const Vehicle_SolverInfo *_solverinfo,
    const Vehicle_Problem *_prob,
    const double _var[...],
    const double _scaling_factor[...],
    double _l[...],
    double _jac_l[...],
    double _hess[...])
{
    if (_l != NULL) {
    /* define objective here */
    }
    
    if (_jac_l != NULL) {
    /* define Jacobian here */
    }
    
    if (_hess != NULL) {
    /* define Hessian increment here */
    }
}

The C/C++ function parameters are defined as follows:

• _solverinfo: The solver info struct, which records the thread_id, stage, and iteration at which this function is called
• _prob: Stores the parameters of the problem
• _var: The optimization variables at this stage
• _scaling_factor: The scaling factor for the objective, which needs to be multiplied onto _l, _jac_l, and _hess_l
• _l: The objective function value
• _jac_l: The Jacobian value of the objective
• _hess: The objective's Hessian is added incrementally (_hess += _hess_l)

info

1. When computing _hess in C/C++ functions, the Hessian of the objective or constraint function is computed incrementally rather than by assignment, i.e.: _hess += _hess_l
2. The calling order of functions in the optimization problem during solving is: $l, l_e, f_s, f_e, f, g$
3. Matrices are stored in column-major order: When a matrix is identified as sparse, non-zero elements are stored sequentially in column-major order; when a matrix is dense, the full matrix is stored. The corresponding C/C++ code must be written according to whether the matrix is sparse or not
4. The Hessian matrix only stores the lower triangle (including the diagonal)

Sparse and Dense Matrices

As mentioned above, matrices can be identified as either sparse or dense. When defining external C/C++ functions, you need to write the corresponding code based on the matrix sparsity. How do you determine whether a matrix is sparse or dense? All of the following conditions must be met:

• Only constraint Jacobian matrices may be identified as sparse; objective Jacobian matrices and all Hessian matrices are always dense
• During modeling, the user specified the matrix sparsity (see the sparsity_... parameters of each external interface, e.g., the sparsity_jacobian parameter in external_general_inequality)
• In the generated template code, if the matrix size is smaller than the dense matrix size, then the matrix is identified as sparse

For example, when defining a 4-dimensional inequality constraint and specifying the sparsity of its Jacobian matrix:

vehicle.py

ineq_spy_jac = \
       [[1, 1, 0, 0, 0, 0, 0],
        [1, 0, 1, 0, 0, 0, 0],
        [1, 0, 0, 1, 0, 0, 0],
        [1, 0, 0, 0, 1, 0, 0]]
ineq_spy_hess = \
       [[0, 0, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 0, 0]]
ineq = external_general_inequality(dim = 4, sparsity_jacobian=ineq_spy_jac, sparsity_hessian=ineq_spy_hess)
prob.inequality(ineq)

And in the generated C/C++ template code, the Jacobian matrix is identified as sparse (the function parameter _jac_g has a size of 8, which is smaller than the dense matrix size of 28 (4-dimensional constraint, 7-dimensional variable)):

vehicle_prob_ext.c.template

void vehicle_g(const Vehicle_SolverInfo *_solverinfo, 
               const Vehicle_Problem *_prob, 
               const double _var[7], 
               const double _zg[4], 
               double _g[4], double _jac_g[8], double _hess[28])
{
    if (_g != NULL) {
    /* define constraint here */
    }
    
    if (_jac_g != NULL) {
    /* define Jacobian here */
    }
    
    if (_hess != NULL) {
    /* define Hessian increment here */
    }
}

warning

Even after specifying matrix sparsity during modeling, the matrix may not necessarily be identified as sparse. You need to check the function parameter sizes in the generated C/C++ template code to determine whether the matrix is identified as sparse.

Passing Custom Data

When a function is an external C/C++ function, custom data can be passed through the _prob->external_data interface. The interface is defined as follows:

typedef struct Vehicle_Problem {
    ...
    void *external_data; /* user defined data */
} Vehicle_Problem;

An example of passing custom data through this interface is as follows:

vehicle_demo.c

Vehicle_Problem prob;
MyData data;
data.a = 1.0;
data.b = 2.0;
prob.external_data = &data;
...

Then in the external C/C++ function, custom data can be accessed through _prob->external_data, for example:

vehicle_userdef_prob_ext.c

void vehicle_l(
    const Vehicle_SolverInfo *_solverinfo,
    const Vehicle_Problem *_prob,
    const double _var[...],
    const double _scaling_factor[...],
    double _out[...],
    double _jac_out[...],
    double _hess[...])
{
    MyData* data = (MyData*)(_prob->external_data);
    double a = data->a;
    double b = data->b;
    ...
}