外部C/C++函数支持

上述章节主要介绍了如何通过符号表达式显式地定义问题, 当优化问题不能通过符号定义时 (例如约束为神经网络或者复杂的机器人动力学模型时) , 可通过外部C/C++函数直接定义.

自定义函数

OPTIMake通过external接口支持以下函数的C/C++定义:

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

在定义外部函数时, 需要提供该函数以及其Jacobian实现. 在某些场景中 (根据option中hessian_approximation的配置), 需要提供Hessian实现:

hessian_approximation配置	外部objective是否提供Hessian	外部约束是否提供Hessian
'gauss-newton'	是	否
'exact'	是	是
'bfgs'	否	否

在使用external接口定义objective或约束为外部函数后, 代码生成功能会自动生成该external函数的C/C++模板.

以下为objective的外部C/C++函数的模板:

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 */
    }
}

其中, C/C++函数入参的定义如下:

• _solverinfo: solver info的结构体, 记录了该函数被调用的thread_id, stage以及迭代
• _prob: 存储了该问题的参数
• _var: 该stage的优化变量
• _scaling_factor: objective的scaling factor, 需要在_l, _jac_l, _hess_l上乘以该factor
• _l: objective的函数值
• _jac_l: objective的Jacobian函数值
• _hess: 将objective的Hessian通过增量方式加入 (_hess += _hess_l)

信息

1. 在计算C/C++函数中的_hess时, objective或约束函数的Hessian通过增量而不是赋值计算, 即: _hess += _hess_l
2. 优化问题中函数在求解过程中的调用顺序为: $l, l_e, f_s, f_e, f, g$
3. 矩阵通过列优先存储: 当矩阵被识别成稀疏矩阵时, 按列优先的方式依次存储非零元素; 当矩阵为稠密矩阵时, 存储完整矩阵. 需要根据矩阵是否稀疏来编写对应的C/C++代码
4. Hessian矩阵只存储下三角 (包含对角线)

稀疏与稠密矩阵

在上面的介绍中, 我们提到矩阵可以被识别成稀疏矩阵或稠密矩阵. 在定义外部C/C++函数时, 需要根据矩阵的稀疏性编写对应的代码. 那么如何判断一个矩阵是稀疏还是稠密矩阵呢? 我们有以下原则:

• 只有约束的Jacobian矩阵可能被识别成稀疏矩阵, objective的Jacobian矩阵和所有Hessian矩阵均为稠密矩阵
• 在建模时, 用户指定了矩阵的稀疏性 (见各external接口的sparsity_...参数, 例如external_general_inequality中的sparsity_jacobian参数)
• 在生成的模板代码中, 矩阵的大小小于稠密矩阵的大小, 则该矩阵被识别成稀疏矩阵

例如, 当定义了一个4维的不等式约束, 并且指定了该不等式约束的Jacobian矩阵的稀疏性:

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)

并且在生成的C/C++模板代码中, 该Jacobian矩阵被识别成稀疏矩阵 (函数参数_jac_g的大小为8, 小于稠密矩阵的大小28 (4维约束, 7维变量)):

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 */
    }
}

传入自定义数据

当一个函数是外部C/C++函数时, 该函数可以通过_prob->external_data接口传递自定义数据. 该接口的定义如下:

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

一个通过该接口传递自定义数据的例子如下:

vehicle_demo.c

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

然后在外部C/C++函数中通过_prob->external_data可获取自定义数据, 例如

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;
    ...
}