调查线性不可行性

此示例说明如何调查导致问题不可行的线性约束。有关这些技术的更多详细信息，请参阅 Chinneck [1] 和 [2]。

如果线性约束导致问题不可行，则您可能需要找到不可行的约束子集，但删除该子集中的任何一个成员都会使该子集的其余部分变得可行。这种子集的名称是不可约不可行约束集，缩写为 IIS。相反，您可能希望找到可行的最大基数约束集。这个子集被称为最大可行子集，缩写为 MaxFS。这两个概念相关，但并不相同。一个问题可能有许多不同的 IIS，其中一些具有不同的基数。

此示例显示了查找 IIS 的两种方法，以及获取可行约束集的两种方法。该示例假设所有给定的边界都是正确的，这意味着 lb 和 ub 参量没有错误。

不可行示例

创建一个随机矩阵 A，表示大小为 150×15 的线性不等式。将相应的向量 b 设置为具有 10 个条目的向量，并将这些值的 5% 更改为 -10。

N = 15;
rng default
A = randn([10*N,N]);
b = 10*ones(size(A,1),1);
Aeq = [];
beq = [];
b(rand(size(b)) <= 0.05) = -10;
f = ones(N,1);
lb = -f;
ub = f;

检查 problem 是否不可行。

[x,fval,exitflag,output,lambda] = linprog(f,A,b,Aeq,beq,lb,ub);

No feasible solution found.

Linprog stopped because no point satisfies the constraints.

删除过滤器

要识别 IIS，请执行以下步骤。给定一组编号为 1 到 N 的线性约束，其中所有问题约束都是不可行的：

对于从 1 到 i 的每个 N：

暂时从问题中删除约束 i。
测试所得到问题的可行性。
如果在没有约束 i 的情况下问题是可行的，则将约束 i 返回到问题。
如果问题在没有约束 i 的情况下可行，则不要将约束 i 返回到问题中

继续下一个 i（直到值 N）。

在此过程结束时，问题中剩余的约束形成了 IIS。

有关实现此过程的 MATLAB® 代码，请参阅此示例末尾的 deletionfilter 辅助函数。

注意：如果您使用此示例的实时脚本文件，则 deletionfilter 函数已包含在文件末尾。否则，您需要在 .m 文件末尾创建此函数或将其作为文件添加到 MATLAB 路径上。本示例后面使用的其他辅助函数也是如此。

查看 deletionfilter 对示例数据的影响。

[ineqs,eqs,ncall] = deletionfilter(A,b,Aeq,beq,lb,ub);

该问题没有等式约束。找到不等式约束的索引和 b(iis) 的值。

iis = find(ineqs)

iis = 
114

b(iis)

ans = 
-10

仅一个不等式约束与边界约束一起就会导致问题不可行。约束为

A(iis,:)*x <= b(iis).

为什么这个约束与边界一起是不可行的？求出 A 该行的绝对值之和。

disp(sum(abs(A(iis,:))))

    8.4864

由于边界，x 向量的值介于 -1 和 1 之间，因此 A(iis,:)*x 不能小于 b(iis) = -10。

linprog 执行了多少次 deletionfilter 调用？

disp(ncall)

该问题有 150 个线性约束，因此该函数调用了 linprog 150 次。

弹性过滤器

作为检查每个约束的删除过滤器的替代，请尝试弹性过滤器。该过滤器的工作原理如下。

首先，允许每个约束 i 被正量 e(i) 违反，其中等式约束具有加法和减法正弹性值。

$\begin{array}{l} A_{i n e q} x \leq b_{i n e q} + e \\ A_{e q} x = b_{e q} + e_{1} - e_{2} \end{array}$

接下来，求解相关的线性规划问题 (LP)

$\min_{x, e} \sum e_{i}$

受到列出的约束并符合 $e_{i} \geq 0$ 。

如果相关 LP 有一个解，则删除所有具有严格正相关的 $e_{i}$ 的约束，并将这些约束记录在索引列表（潜在 IIS 成员）中。返回上一步来求解新的、简化的相关 LP。
如果相关的 LP 没有解（不可行）或者没有严格正相关的 $e_{i}$ ，则退出过滤器。

弹性过滤器可以比删除过滤器在更少的迭代中退出，因为它可以一次将许多索引带入 IIS，并且可以在不遍历整个索引列表的情况下停止。但是，该问题比原始问题具有更多的变量，并且其生成的索引列表可能比 IIS 更大。要在运行弹性过滤器后找到 IIS，请对结果运行删除过滤器。

有关实现此过滤器的 MATLAB® 代码，请参阅此示例末尾的 elasticfilter 辅助函数。

查看 elasticfilter 对示例数据的影响。

[ineqselastic,eqselastic,ncall] = ...
    elasticfilter(A,b,Aeq,beq,lb,ub);

该问题没有等式约束。找到不等式约束的索引。

iiselastic = find(ineqselastic)

iiselastic = 5×1

     2
    60
    82
    97
   114

弹性 IIS 列出了五个约束，而删除过滤器只发现了一个。对返回的集合运行删除过滤器以查找真正的 IIS。

A_1 = A(ineqselastic > 0,:);
b_1 = b(ineqselastic > 0);
[dineq_iis,deq_iis,ncall2] = deletionfilter(A_1,b_1,Aeq,beq,lb,ub);
iiselasticdeletion = find(dineq_iis)

iiselasticdeletion = 
5

弹性滤波结果中的第五个约束，不等式 114，是 IIS。该结果与删除过滤器的答案一致。这两种方法的区别在于，组合弹性和删除过滤器方法使用的 linprog 调用要少得多。显示弹性过滤器和其后的删除过滤器使用的 linprog 调用总数。

disp(ncall + ncall2)

循环删除 IIS

通常，获得单个 IIS 并不能使您找到优化问题失败的所有原因。要纠正一个不可行问题，您可以反复查找 IIS 并将其从问题中删除，直到问题变得可行。

下面的代码显示了如何从问题中一次删除一个 IIS，直到问题变得可行。在算法删除任何约束之前，代码使用索引技术来跟踪约束在原始问题中的位置。

代码使用布尔向量 activeA 来表示 A 矩阵的当前约束（行），使用布尔向量 activeAeq 来表示 Aeq 矩阵的当前约束，从而跟踪问题中的原始变量。当添加或删除约束时，代码会索引到 A 或 Aeq，这样即使约束数量发生变化，行号也不会改变。

运行此代码将返回 idx2，即 activeA 中非零元素索引的向量：

idx2 = find(activeA)

假设 var 是一个布尔向量，其长度与 idx2 相同。然后，

idx2(find(var))

将 var 表示为原始问题变量的索引。通过这种方式，索引可以采用约束集的子集，只使用较小的子集，并且仍然明确地引用原始问题变量。

opts = optimoptions('linprog','Display',"none");
activeA = true(size(b));
activeAeq = true(size(beq));
[~,~,exitflag] = linprog(f,A,b,Aeq,beq,lb,ub,opts);
ncl = 1;
while exitflag < 0
    [ineqselastic,eqselastic,ncall] = ...
        elasticfilter(A(activeA,:),b(activeA),Aeq(activeAeq,:),beq(activeAeq),lb,ub);
    ncl = ncl + ncall;
    idxaa = find(activeA);
    idxae = find(activeAeq);
    tmpa = idxaa(find(ineqselastic));
    tmpae = idxae(find(eqselastic));
    AA = A(tmpa,:);
    bb = b(tmpa);
    AE = Aeq(tmpae,:);
    be = beq(tmpae);
    [ineqs,eqs,ncall] = ...
        deletionfilter(AA,bb,AE,be,lb,ub);
    ncl = ncl + ncall;
    activeA(tmpa(ineqs)) = false;
    activeAeq(tmpae(eqs)) = false;
    disp(['Removed inequalities ',int2str((tmpa(ineqs))'),' and equalities ',int2str((tmpae(eqs))')])
    [~,~,exitflag] = ...
        linprog(f,A(activeA,:),b(activeA),Aeq(activeAeq,:),beq(activeAeq),lb,ub,opts);
    ncl = ncl + 1;
end

Removed inequalities 114 and equalities 
Removed inequalities 97 and equalities 
Removed inequalities 64  82 and equalities 
Removed inequalities 60 and equalities

fprintf('Number of linprog calls: %d\n',ncl)

Number of linprog calls: 28

请注意，循环同时消除了不等式 64 和 82，这表明这两个约束形成了 IIS。

查找 MaxFS

获得可行约束集的另一种方法是直接找到 MaxFS。正如 Chinneck [1] 所解释的，寻找 MaxFS 是一个 NP 完全问题，这意味着该问题不一定存在寻找 MaxFS 的有效算法。然而，Chinneck 提出了一些可以有效工作的算法。

使用 Chinneck 算法 7.3 来查找约束的覆盖集，当将其删除时，可得到可行集。该算法在本示例末尾的 generatecover 辅助函数中实现。

[coversetineq,coverseteq,nlp] = generatecover(A,b,Aeq,beq,lb,ub)

coversetineq = 5×1

   114
    97
    60
    82
     2

coverseteq =

     []

nlp = 
40

消除这些约束并求解 LP。

usemeineq = true(size(b));
usemeineq(coversetineq) = false; % Remove inequality constraints
usemeeq = true(size(beq));
usemeeq(coverseteq) = false; % Remove equality constraints
[xs,fvals,exitflags] = ...
    linprog(f,A(usemeineq,:),b(usemeineq),Aeq(usemeeq),beq(usemeeq),lb,ub);

Optimal solution found.

请注意，封面集与 Elastic Filter 中的 iiselastic 集完全相同。一般来说，弹性过滤器会发现过大的覆盖集。Chinneck 算法 7.3 从弹性滤波结果开始，然后仅保留必要的约束。

Chinneck 算法 7.3 需要调用 40 次 linprog 才能完成 MaxFS 的计算。这个数字比之前在循环删除 IIS 的过程中使用的 28 次调用要多一点。

另外，请注意循环中消除的不等式与算法 7.3 消除的不等式并不完全相同。循环消除了不等式 114、97、82、60 和 64，而算法 7.3 消除了不等式 114、97、82、60 和 2。检查不等式 82 和 64 是否形成 IIS（如循环中删除 IIS 中所示），以及不等式 82 和 2 是否也形成 IIS。

usemeineq = false(size(b));
usemeineq([82,64]) = true;
ineqs = deletionfilter(A(usemeineq,:),b(usemeineq),Aeq,beq,lb,ub);
disp(ineqs)

   1
   1

usemeineq = false(size(b));
usemeineq([82,2]) = true;
ineqs = deletionfilter(A(usemeineq,:),b(usemeineq),Aeq,beq,lb,ub);
disp(ineqs)

   1
   1

参考资料

[1] Chinneck, J. W. Feasibility and Infeasibility in Optimization:Algorithms and Computational Methods.Springer, 2008.

[2] Chinneck, J. W."Feasibility and Infeasibility in Optimization."Tutorial for CP-AI-OR-07, Brussels, Belgium.Available at https://www.sce.carleton.ca/faculty/chinneck/docs/CPAIOR07InfeasibilityTutorial.pdf.

辅助函数

以下代码创建 deletionfilter 辅助函数。

function [ineq_iis,eq_iis,ncalls] = deletionfilter(Aineq,bineq,Aeq,beq,lb,ub)
ncalls = 0;
[mi,n] = size(Aineq); % Number of variables and linear inequality constraints
f = zeros(1,n);
me = size(Aeq,1); % Number of linear equality constraints
opts = optimoptions("linprog","Algorithm","dual-simplex","Display","none");

ineq_iis = true(mi,1); % Start with all inequalities in the problem
eq_iis = true(me,1); % Start with all equalities in the problem

for i=1:mi
    ineq_iis(i) = 0; % Remove inequality i
    [~,~,exitflag] = linprog(f,Aineq(ineq_iis,:),bineq(ineq_iis),...
                                Aeq,beq,lb,ub,opts);
    ncalls = ncalls + 1;
    if exitflag == 1 % If now feasible
        ineq_iis(i) = 1; % Return i to the problem
    end
end
for i=1:me
    eq_iis(i) = 0; % Remove equality i
    [~,~,exitflag] = linprog(f,Aineq,bineq,...
                                Aeq(eq_iis,:),beq(eq_iis),lb,ub,opts);
    ncalls = ncalls + 1;
    if exitflag == 1 % If now feasible
        eq_iis(i) = 1; % Return i to the problem
    end
end
end

以下代码创建 elasticfilter 辅助函数。

function [ineq_iis,eq_iis,ncalls,fval0] = elasticfilter(Aineq,bineq,Aeq,beq,lb,ub)
ncalls = 0;
[mi,n] = size(Aineq); % Number of variables and linear inequality constraints
me = size(Aeq,1);
Aineq_r = [Aineq -1.0*eye(mi) zeros(mi,2*me)];
Aeq_r = [Aeq zeros(me,mi) eye(me) -1.0*eye(me)]; % Two slacks for each equality constraint
lb_r = [lb(:); zeros(mi+2*me,1)];
ub_r = [ub(:); inf(mi+2*me,1)];
ineq_slack_offset = n;
eq_pos_slack_offset = n + mi;
eq_neg_slack_offset = n + mi + me;
f = [zeros(1,n) ones(1,mi+2*me)];
opts = optimoptions("linprog","Algorithm","dual-simplex","Display","none");
tol = 1e-10;

ineq_iis = false(mi,1);
eq_iis = false(me,1);
[x,fval,exitflag] = linprog(f,Aineq_r,bineq,Aeq_r,beq,lb_r,ub_r,opts);
fval0 = fval;
ncalls = ncalls + 1;
while exitflag == 1 && fval > tol % Feasible and some slacks are nonzero
    c = 0;
    for i = 1:mi        
        j = ineq_slack_offset+i;
        if x(j) > tol
            ub_r(j) = 0.0;
            ineq_iis(i) = true;
            c = c+1;
        end
    end
    for i = 1:me
        j = eq_pos_slack_offset+i;
        if x(j) > tol            
            ub_r(j) = 0.0;
            eq_iis(i) = true;
            c = c+1;
        end
    end
    for i = 1:me
        j = eq_neg_slack_offset+i;
        if x(j) > tol            
            ub_r(j) = 0.0;
            eq_iis(i) = true;
            c = c+1;
        end
    end
    [x,fval,exitflag] = linprog(f,Aineq_r,bineq,Aeq_r,beq,lb_r,ub_r,opts);
    if fval > 0
        fval0 = fval;
    end
    ncalls = ncalls + 1;
end 
end

以下代码创建 generatecover 辅助函数。该代码使用与循环中删除 IIS 代码相同的索引技术来跟踪约束。

function [coversetineq,coverseteq,nlp] = generatecover(Aineq,bineq,Aeq,beq,lb,ub)
% Returns the cover set of linear inequalities, the cover set of linear
% equalities, and the total number of calls to linprog.
% Adapted from Chinneck [1] Algorithm 7.3. Step numbers are from this book.
coversetineq = [];
coverseteq = [];
activeA = true(size(bineq));
activeAeq = true(size(beq));
% Step 1 of Algorithm 7.3
[ineq_iis,eq_iis,ncalls] = elasticfilter(Aineq,bineq,Aeq,beq,lb,ub);
nlp = ncalls;
ninf = sum(ineq_iis(:)) + sum(eq_iis(:));
if ninf == 1
    coversetineq = ineq_iis;
    coverseteq = eq_iis;
    return
end
holdsetineq = find(ineq_iis);
holdseteq = find(eq_iis);
candidateineq = holdsetineq;
candidateeq = holdseteq;
% Step 2 of Algorithm 7.3
while sum(candidateineq(:)) + sum(candidateeq(:)) > 0
    minsinf = inf;
    ineqflag = 0;
    for i = 1:length(candidateineq(:))
        activeA(candidateineq(i)) = false;
        idx2 = find(activeA);
        idx2eq = find(activeAeq);
        [ineq_iis,eq_iis,ncalls,fval] = elasticfilter(Aineq(activeA,:),bineq(activeA),Aeq(activeAeq,:),beq(activeAeq),lb,ub);
        nlp = nlp + ncalls;
        ineq_iis = idx2(find(ineq_iis));
        eq_iis = idx2eq(find(eq_iis));
        if fval == 0
            coversetineq = [coversetineq;candidateineq(i)];
            return
        end
        if fval < minsinf
            ineqflag = 1;
            winner = candidateineq(i);
            minsinf = fval;
            holdsetineq = ineq_iis;
            if numel(ineq_iis(:)) + numel(eq_iis(:)) == 1
                nextwinner = ineq_iis;
                nextwinner2 = eq_iis;
                nextwinner = [nextwinnner,nextwinner2];
            else
                nextwinner = [];
            end
        end
        activeA(candidateineq(i)) = true;
    end
    for i = 1:length(candidateeq(:))
        activeAeq(candidateeq(i)) = false;
        idx2 = find(activeA);
        idx2eq = find(activeAeq);
        [ineq_iis,eq_iis,ncalls,fval] = elasticfilter(Aineq(activeA,:),bineq(activeA),Aeq(activeAeq,:),beq(activeAeq),lb,ub);
        nlp = nlp + ncalls;
        ineq_iis = idx2(find(ineq_iis));
        eq_iis = idx2eq(find(eq_iis));
        if fval == 0
            coverseteq = [coverseteq;candidateeq(i)];
            return
        end
        if fval < minsinf
            ineqflag = -1;
            winner = candidateeq(i);
            minsinf = fval;
            holdseteq = eq_iis;
            if numel(ineq_iis(:)) + numel(eq_iis(:)) == 1
                nextwinner = ineq_iis;
                nextwinner2 = eq_iis;
                nextwinner = [nextwinnner,nextwinner2];
            else
                nextwinner = [];
            end
        end
        activeAeq(candidateeq(i)) = true;
    end
% Step 3 of Algorithm 7.3
    if ineqflag == 1
        coversetineq = [coversetineq;winner];
        activeA(winner) = false;
        if nextwinner
            coversetineq = [coversetineq;nextwinner];
            return
        end
    end
    if ineqflag == -1
        coverseteq = [coverseteq;winner];
        activeAeq(winner) = false;
        if nextwinner
            coverseteq = [coverseteq;nextwinner];
            return
        end
    end
    candidateineq = holdsetineq;
    candidateeq = holdseteq;
end
end

另请参阅

linprog