如何利用DecomposeLumpedParameters重构MathematicalProgram决策变量？

核心思路 #

借助Drake的符号表达式替换工具，结合DecomposeLumpedParameters自动识别问题中的集总乘积项，将这些项替换为新的决策变量，重构优化问题以实现新变量下的线性化，最终直接求解并提取集总参数的值。

具体实现步骤 #

• 提取集总项：收集原问题中所有目标函数和约束的表达式，传入DecomposeLumpedParameters函数，得到所有集总乘积项与对应变量组合的映射关系。
• 定义新决策变量：为每个识别出的集总项创建新的连续决策变量，建立集总表达式到新变量的映射。
• 重构优化问题：遍历原问题的目标和约束，使用Substitute函数将所有集总项替换为对应的新变量；若需严格保证原变量与新变量的等价性，可添加m1 == a*b这类约束（但会引入非线性，需根据需求选择）。
• 求解并提取结果：调用求解器求解重构后的问题，直接获取新变量（即集总参数）的最优值。

代码示例 #

#include "drake/solvers/mathematical_program.h"
#include "drake/symbolic/symbolic.h"
#include "fmt/core.h"

using namespace drake;
using namespace drake::solvers;
using namespace drake::symbolic;

int main() {
    // 构建原优化问题
    MathematicalProgram prog;
    const Variable a = prog.NewContinuousVariable("a");
    const Variable b = prog.NewContinuousVariable("b");
    const Variable c = prog.NewContinuousVariable("c");

    // 添加仅含集总项的目标与约束
    prog.AddCost(2*a*b + 3*b*c);
    prog.AddConstraint(a*b + b*c <= 10);
    prog.AddConstraint(a*b >= 2);

    // 收集所有待分析的表达式
    std::vector<Expression> all_exprs;
    for (const auto& cost : prog.GetAllCosts()) {
        all_exprs.push_back(cost->evaluator()->Eval({}));
    }
    for (const auto& constraint : prog.GetAllConstraints()) {
        for (int i = 0; i < constraint->num_constraints(); ++i) {
            all_exprs.push_back(constraint->evaluator()->Eval({})[i]);
        }
    }

    // 自动识别集总项
    const auto lumped_map = DecomposeLumpedParameters(all_exprs);

    // 创建新变量并建立映射
    std::map<Expression, Variable> lumped_to_new_var;
    int var_idx = 1;
    for (const auto& [lumped_expr, _] : lumped_map) {
        const Variable new_var = prog.NewContinuousVariable(fmt::format("m{}", var_idx++));
        lumped_to_new_var[lumped_expr] = new_var;
    }

    // 替换目标函数
    prog.ClearCosts();
    const Expression original_cost = 2*a*b + 3*b*c;
    prog.AddCost(Substitute(original_cost, lumped_to_new_var));

    // 替换约束
    prog.ClearConstraints();
    const Expression constr1 = a*b + b*c;
    prog.AddConstraint(Substitute(constr1, lumped_to_new_var) <= 10);
    const Expression constr2 = a*b;
    prog.AddConstraint(Substitute(constr2, lumped_to_new_var) >= 2);

    // 求解问题
    const auto result = Solve(prog);
    if (result.is_success()) {
        std::cout << "m1 (ab) = " << result.GetSolution(lumped_to_new_var.at(a*b)) << "\n";
        std::cout << "m2 (bc) = " << result.GetSolution(lumped_to_new_var.at(b*c)) << "\n";
    }

    return 0;
}

关键注意事项 #

• 若原问题存在非集总形式的非线性项，替换后无法保证问题线性，需提前确认所有非线性仅来自集总乘积项。
• 添加m == a*b这类等价约束会让问题重回非线性，若目标是严格线性化，需根据问题场景选择是否跳过这类约束（比如问题本身仅需集总参数满足线性约束，无需关联原变量）。
• DecomposeLumpedParameters支持识别任意阶数的变量乘积项，不限于二元乘积。

内容的提问来源于stack exchange，提问作者Nicholas Pfaff