带整数约束的曲线线性组合拟合及CPLEX求解报错咨询

问题本质梳理 #

你遇到的是**混合整数二次规划（MIQP）**问题：最小二乘的目标函数是二次的，而你需要的「限制拟合曲线数量」「要求所选曲线连续相邻」这类约束，都需要引入0-1整数变量来实现——普通的lsqlin只能处理连续变量的约束（比如你已经实现的非负限制），没法直接搞定整数约束，所以得换用支持MIQP的求解器，比如CPLEX、Gurobi，或者开源的SCIP。

先解决你当前CPLEX示例的报错 #

你用cplexlsqnonneglin时碰到的CPLEX Error 5002: %s is not convex，原因很明确：这个函数专门处理非负连续变量的凸最小二乘问题，要求C'*C是半正定矩阵。如果你的C存在线性相关的列（列不满秩），就会导致矩阵不满足凸性要求。

不过更重要的是：cplexlsqnonneglin根本不支持整数约束！如果你的最终目标是加整数约束，别在这个函数上浪费时间，直接把问题建模成MIQP才是正确方向。

正确建模思路（带整数约束） #

假设我们要拟合d ≈ C*x，其中x是曲线的系数向量。要实现你要的约束，需要引入0-1整数变量z_j：z_j=1表示选中第j条曲线，z_j=0表示不选。

1. 限制拟合曲线数量 #

比如最多选k条曲线，需要添加两个约束：

• x_j ≤ M * z_j：M是一个足够大的正数，确保当z_j=0时，x_j被强制为0（相当于不选这条曲线）
• sum(z_j) ≤ k：限制选中的曲线总数不超过k
• 保留你之前的非负约束：x_j ≥ 0

2. 要求所选曲线连续相邻 #

如果曲线是按顺序排列的，连续相邻的约束可以转化为线性约束：对于所有j，z_j - z_{j+1} + z_{j+2} ≤ 1。这个约束的逻辑是：如果第j和j+2条曲线都被选中，那么中间的j+1条必须也被选中，避免出现“跳选”的情况。

最终的目标函数是最小化0.5*(C*x - d)'*(C*x - d)，这是一个二次目标，加上上面的线性约束和整数变量，就构成了标准的MIQP问题，可以用CPLEX的cplexmiqp函数求解。

针对你的场景的代码示例 #

下面是一个可参考的MIQP建模代码框架，适配你的需求：

function [] = miqp_curve_fit_example()
    load('C_n_d_2.mat');
    n = size(C,2); % 曲线的总数量
    max_selected = 5; % 你想要限制的最多选中曲线数，按需调整
    % 设定M值：要足够大，确保能约束x_j为0，但也别太大避免数值不稳定
    M = max(max(abs(C))) * norm(d); 
    
    % 构建MIQP的目标函数：0.5*x'*Q*x + f'*x
    Q = C' * C;
    f = -C' * d;
    
    % 构建约束条件 Ax ≤ b
    % 约束1：x_j ≤ M*z_j → x_j - M*z_j ≤ 0
    A1 = [eye(n), -M*eye(n)];
    b1 = zeros(n, 1);
    % 约束2：选中的曲线总数不超过max_selected
    A2 = [zeros(1, n), ones(1, n)];
    b2 = max_selected;
    % 如果需要连续相邻约束，添加下面的代码块
    % A3 = sparse(n-2, 2*n);
    % for j = 1:n-2
    %     A3(j, j) = 1;
    %     A3(j, j+1) = -1;
    %     A3(j, j+2) = 1;
    %     A3(j, n+j) = -1;
    %     A3(j, n+j+1) = 1;
    %     A3(j, n+j+2) = -1;
    % end
    % b3 = ones(n-2, 1);
    % 合并约束（如果加连续约束，修改下面两行的参数）
    A = [A1; A2];
    b = [b1; b2];
    
    % 变量上下界：x≥0，z是0-1变量
    lb = [zeros(n, 1); zeros(n, 1)];
    ub = [inf(n, 1); ones(n, 1)];
    
    % 变量类型：前n个是连续变量（'C'），后n个是0-1整数变量（'B'）
    ctype = [repmat('C', 1, n), repmat('B', 1, n)];
    
    % CPLEX求解设置
    options = cplexoptimset('Display', 'on');
    [x_z, fval, exitflag, output] = cplexmiqp(Q, f, A, b, [], [], lb, ub, ctype, [], options);
    
    % 提取结果
    x_opt = x_z(1:n); % 拟合系数
    z_opt = x_z(n+1:end); % 选中曲线的0-1标记
    disp('选中的曲线索引：');
    disp(find(z_opt == 1));
    disp('拟合得到的系数：');
    disp(x_opt);
end

额外小贴士 #

• 如果你用开源工具，可以试试SCIP，它对MIQP的支持也很好，不需要商业授权。
• M值的选择很关键：可以用max(max(abs(C))) * norm(d)来估算，既能保证约束有效，又不会引入数值问题。
• 连续相邻约束的逻辑可以根据你的具体需求调整，比如如果允许多个连续区间，上面的约束也能适用。

内容的提问来源于stack exchange，提问作者Jesse RJ