版本：v0.7.0 beta

可微参数

OPTIMake支持计算最优解关于parameter的梯度 (参数灵敏度). 在建模时将parameter声明为可微参数后, 生成的代码中会额外提供backward接口, 可用于:

计算MPC的局部反馈增益, 例如首个控制量关于初始状态的导数 $\partial u_0 / \partial x_0$
分析最优解关于参考轨迹, 权重等参数的灵敏度
将求解器作为可微的优化层, 通过梯度下降学习模型中的参数 (例如目标函数权重, 约束参数)

计算原理

记优化变量为 $v$ , 可微参数为 $\theta$ , 最优解为 $v^*(\theta)$ . 对于用户关心的标量函数 $c(v)$ , 用户提供其关于最优解的梯度 $\nabla_v c$ , backward接口计算 $c$ 关于可微参数的梯度 (向量-雅可比积):

\nabla_\theta c = \left(\frac{\partial v^*}{\partial \theta}\right)^T \nabla_v c

backward基于最优解处KKT条件的隐式微分, 不对求解过程中的迭代, 线搜索等步骤进行反向传播, 每次调用的计算量约为一次KKT系统的求解.

信息

若 $c$ 还直接依赖于参数 $\theta$ (即 $c(v, \theta)$ ), 用户需要自行加上 $v$ 固定时 $c$ 关于 $\theta$ 的偏导数.

建模

在parameter或parameters接口中设置differentiable=True, 即可将parameter声明为可微参数. stage-independent与stage-dependent parameter均支持:

from optimake import *

prob = multi_stage_problem('tracking', 20)

p = prob.variable('p')
v = prob.variable('v')
a = prob.variable('a')

# 可微参数
p0, v0 = prob.parameters(['p0', 'v0'], differentiable=True)
weight = prob.parameter('weight', differentiable=True)
p_ref = prob.parameter('p_ref', stage_dependent=True, differentiable=True)

prob.objective(general_objective(weight * (p - p_ref)**2 + 0.1 * a**2))
prob.equality(differential_equation(state=[p, v], state_dot=[v, a],
                                    stepsize=0.1, discretization_method='erk4'))
prob.start_equality(general_equality([p - p0, v - v0]))

option = codegen_option()
option.default_tolerance_level = 'high'
code_generator().codegen(prob, option)

differentiable默认为False. 可微参数仍然是普通的parameter, 在param与param_stage中的设置方式与index均不变
只要模型中存在可微参数, 代码生成时即会生成backward接口, 无需额外的代码生成选项
求解算法'pdipm'与'sqp'均支持

生成的接口

以上述tracking模型为例, 与可微参数相关的生成内容如下:

名称	说明
`problem.grad_var[N][DIM_VAR]`	输入, $\nabla_v c$ , 变量顺序与`output.primal.var`一致, 通过`TRACKING_VAR_*`索引
`output.grad_param[DIM_DP]`	输出, $c$ 关于stage-independent可微参数的梯度
`output.grad_param_stage[N][DIM_DPS]`	输出, $c$ 关于各stage的stage-dependent可微参数的梯度
`TRACKING_DIM_DP`, `TRACKING_DIM_DPS`	stage-independent与stage-dependent可微参数的数量
`TRACKING_DIFF_<NAME>`	stage-independent可微参数在`grad_param`中的index, 例如`TRACKING_DIFF_WEIGHT`
`TRACKING_DIFF_STAGE_<NAME>`	stage-dependent可微参数在`grad_param_stage[i]`中的index, 例如`TRACKING_DIFF_STAGE_P_REF`

注意

梯度数组只包含可微参数, 其index与param/param_stage中的index不同, 读取梯度时请使用DIFF宏 (例如TRACKING_DIFF_WEIGHT), 而不是PARAM宏 (例如TRACKING_PARAM_WEIGHT). 当某一类可微参数的数量为0时, 不会生成对应的grad_param或grad_param_stage.

backward接口定义如下:

int tracking_backward(const Tracking_Problem *const prob,
                      const Tracking_Option *const option,
                      Tracking_WorkSpace *const workspace,
                      Tracking_Output *const output, int refactor);

prob, option, workspace: 与前一次solve使用的相同, 其中prob->grad_var为输入的 $\nabla_v c$
output: 写入grad_param与grad_param_stage, 不修改output中的primal变量等求解结果
refactor:
- 1: 在最终解处重新计算导数并重新分解KKT系统, 再计算梯度
- 0: 复用已有的KKT分解 (来自前一次求解或最近一次refactor=1的调用), 只求解一次线性系统
返回值: 1表示成功; 0表示失败, 此时梯度数组被清零, 并打印失败原因

使用方法

调用流程为: solve → 设置grad_var → backward → 读取梯度. 下面的例子分别计算首个控制量 $a_0$ 与最后一个stage的位置 $p_{N-1}$ 关于可微参数的梯度:

int solve_status = tracking_solve(&problem, &option, &workspace, &output);
if (solve_status == 1) {
    /* c = a_0 */
    memset(problem.grad_var, 0, sizeof(problem.grad_var));
    problem.grad_var[0][TRACKING_VAR_A] = 1.0;
    /* 每次求解后的第一次调用使用 refactor = 1 */
    if (tracking_backward(&problem, &option, &workspace, &output, 1) == 1) {
        double da0_dp0 = output.grad_param[TRACKING_DIFF_P0];
        double da0_dweight = output.grad_param[TRACKING_DIFF_WEIGHT];
        double da0_dpref_5 = output.grad_param_stage[5][TRACKING_DIFF_STAGE_P_REF];
    }

    /* c = p_{N-1}, 更换 grad_var 后使用 refactor = 0 复用分解 */
    memset(problem.grad_var, 0, sizeof(problem.grad_var));
    problem.grad_var[problem.valid_dim_N - 1][TRACKING_VAR_P] = 1.0;
    if (tracking_backward(&problem, &option, &workspace, &output, 0) == 1) {
        double dpN_dp0 = output.grad_param[TRACKING_DIFF_P0];
    }
}

信息

每次调用backward得到一个标量函数的梯度, 即雅可比矩阵的一行. 若需要完整的雅可比矩阵 (例如 $\partial u_0 / \partial x_0$ ), 可依次将grad_var设置为各分量的单位向量并多次调用backward: 第一次使用refactor=1, 其余使用refactor=0.

注意事项

求解状态: 梯度在收敛的最优解处才有意义, 建议仅在solve_status为1时使用backward的结果. backward本身不检查求解状态, 对于其他求解状态, 返回值仅为当前局部系统的灵敏度估计
调用顺序: 在solve与backward之间, 除grad_var外, 不要修改problem (包括parameter与valid_dim_N), option与workspace. parameter改变后需要重新求解再调用backward
输出: 每次调用backward都会覆盖grad_param与grad_param_stage, solve不会更新这两个数组. valid_dim_N之后的stage以及未在模型中使用的可微参数, 其梯度为0
refactor的选择: 求解结束时保留的KKT分解通常来自最后一次迭代更新之前, 直接使用refactor=0可能引入额外误差. 建议每次求解后的第一次调用使用refactor=1, 之后对同一次求解结果更换grad_var时使用refactor=0
workspace: backward会使用workspace作为计算空间, 但不会修改其中的primal/dual变量, 因此不影响下一次求解的warm start
精度: 梯度精度受求解精度影响, 建议使用较高的求解精度 (例如代码生成选项default_tolerance_level = 'high'). 使用'gauss-newton'或'bfgs'作为Hessian近似时, 梯度基于相应的近似Hessian计算, 为近似梯度; 对精度要求较高时建议使用'exact'
不可微点: 不等式约束在活跃与非活跃之间切换的位置, 以及软约束硬化前后的切换位置, 导数可能不唯一. 当变量被活跃的上下界限制时, 其关于参数的导数为0. 无效化的约束对梯度无贡献
不支持的情况: 以下情况在代码生成时会报错
- 使用外部C/C++函数定义的问题
- 使用'irk4'积分方法的differential_equation
- 占据网格的resolution依赖于可微参数

计算原理​

建模​

生成的接口​

使用方法​

注意事项​

计算原理

建模

生成的接口

使用方法

注意事项