研究线性不可行性

打开实时脚本

这个例子展示了如何研究导致问题不可行的线性约束。有关这些技术的更多细节，请参见Chinneck[1]而且[２]．

如果线性约束导致问题不可行，您可能希望找到不可行约束的子集，但删除子集中的任何成员将使子集的其余部分可行。这种子集的名称是约束的不可约不可行子集，缩写为IIS。相反，您可能希望找到可行的约束的最大基数子集。这个子集叫做a最大可行子集，缩写maxf。这两个概念是相关的，但不是完全相同的。一个问题可以有许多不同的IISs，其中一些具有不同的基数。

此示例展示了查找IIS的两种方法和获取可行约束集的两种方法。该示例假设所有给定的边界都是正确的，即磅而且乌兰巴托参数没有错误。

不可行的例子

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

N = 15;rng默认的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;

检查问题是不可行的。

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

没有可行的解决方案。Linprog停止了，因为没有一个点满足约束条件。

删除过滤器

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

为每一个我从1到N：

暂时取消约束我从问题中。
测试得到的问题的可行性。
如果问题是可行的，没有约束我，返回约束我解决问题。
如果没有约束，问题是不可行的我，不返回约束我解决这个问题

继续下一个我(视价值而定)N）.

在这个过程的最后，保留在问题中的约束形成一个IIS。

关于实现此过程的MATLAB®代码，请参见deletionfilter的辅助功能此示例结束．

请注意:如果在此示例中使用实时脚本文件，则deletionfilter函数已经包含在文件的末尾。否则，您需要在.m文件的末尾创建这个函数，或者将它作为一个文件添加到MATLAB路径中。本例后面使用的其他helper函数也是如此。

看到的效果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)．

为什么这个约束条件和边界是不可行的?求这一行的绝对值的和一个．

disp(总和(abs ((iis,:))))

8.4864

由于有边界xVector的值在-1到1之间，以此类推:一个(iis) * x不能少于b (iis)= -10。

有多少linprog调用了deletionfilter执行?

disp (ncall)

这个问题有150个线性约束条件，所以函数叫做linprog150次。

弹性过滤器

删除筛选器检查每个约束，作为它的替代选择，可以尝试弹性筛选器。这个过滤器的工作原理如下。

首先，允许每个约束我被侵犯:被正数侵犯e(我)，其中等式约束同时具有加式和减式正弹性值。

$\begin{array}{l} {一个}_{我 n e 问} x \leq b_{我 n e 问} + e \\ {一个}_{e 问} x ＝ b_{e 问} + e_{1} - e_{2} \end{array}$

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

$\underset{x ， e}{最小值} \sum e_{我}$

根据所列的约束条件 $e_{我} \geq 0$ ．

如果相关的LP有一个解，则删除所有与严格正相关的约束 $e_{我}$ ，并在索引(潜在的IIS成员)列表中记录这些约束。返回到上一步以解决新的、减少的相关LP。
如果相关的LP没有解(不可行)或没有严格的正相关 $e_{我}$ ，退出过滤。

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

有关实现该过滤器的MATLAB®代码，请参阅elasticfilter的辅助功能此示例结束．

看到的效果elasticfilter在示例数据上。

[ineqselastic, eqselastic ncall] =.．.Aeq elasticfilter (A, b,说真的,磅,乌兰巴托);

这个问题没有等式约束。找到不等式约束的指标。

iselastic = find(非弹性的)

iiselastic =5×12 60 82 97 114

弹性IIS列出了五个约束，而删除筛选器只找到了一个。在返回的集合上运行删除过滤器来查找真正的IIS。

A_1 = A(ineqselastic > 0，:);B_1 = b(非均匀> 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控件的当前约束(行)一个矩阵和布尔向量activeAeq的当前约束条件Aeq矩阵。当添加或删除约束时，代码索引到一个或Aeq这样行号就不会改变，即使约束条件的数量改变了。

运行此代码将返回idx2中的非零元素指标的向量activeA：

idx2 = find(activeA)

假设var是一个布尔向量，长度与idx2．然后

idx2(找到(var))

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

Opts = optimoptions(“linprog”，“显示”，“没有”）;activeA = true(size(b));activeAeq = true(size(beq));[~, ~, exitflag] = linprog (f, A、b Aeq,说真的,磅,乌兰巴托,选择);NCL = 1;而Exitflag < 0 [ineqselastic,eqselastic,ncall] =.．.elasticfilter ((activeA:), b (activeA) Aeq (activeAeq:),说真的(activeAeq)磅,乌兰巴托);NCL = NCL + ncall;idxaa = find(activeA);idxae = find(activeAeq);Tmpa = idxaa(find(ineqselastic));tpae = idxae(find(eqselastic));AA = A(tmpa，:);Bb = b(tmpa);AE = Aeq(tpae，:);Be = beq(tpae); [ineqs,eqs,ncall] =.．.deletionfilter (AA、bb AE,磅,乌兰巴托);NCL = NCL + ncall;activeA(tmpa(ineqs)) = false;activeAeq(tpae (eqs)) = false;disp ([“消除不平等”, int2str ((tmpa (ineqs))”),“和等式”，int2str((tpae (eq))')]) [~，~，exitflag] =.．.linprog (f (activeA:)、b (activeA) Aeq (activeAeq:),说真的(activeAeq)磅,乌兰巴托,选择);NCL = NCL + 1;结束

消除不平等114与平等消除不平等97与平等消除不平等64 82与平等消除不平等60与平等

流(' linprog调用数:%d\n'ncl)

linprog调用数:28

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

找到maxf

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

使用Chinneck的7.3算法找到一个覆盖集当去掉约束时，给出一个可行集。算法实现在generatecover的辅助功能此示例结束．

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

coversetineq =5×1114 97 60 82

Coverseteq = []

NLP = 40

消除这些约束，解决LP。

Usemeineq = true(size(b));Usemeineq (coversetineq) = false;消除不平等约束Usemeeq = true(size(beq));Usemeeq (coverseteq) = false;移除等式约束[x, fvals exitflags] =.．.linprog (f (usemeineq:)、b (usemeineq) Aeq (usemeeq),说真的(usemeeq)磅,乌兰巴托);

找到最优解。

注意，封面集是完全相同的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)

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

参考文献

[1] Chinneck, j.w.。优化的可行性与不可行性:算法与计算方法。施普林格,2008年。

[2] Chinneck, j.w.。"优化的可行性和不可行性"CP-AI-OR-07教程，布鲁塞尔，比利时。可以在https://www.sce.carleton.ca/faculty/chinneck/docs/CPAIOR07InfeasibilityTutorial.pdf．

辅助函数

此代码创建deletionfilterhelper函数。

函数[ineq_iis,eq_iis,ncalls] = deletionfilter(Aineq,bineq,Aeq,beq,lb,ub) ncalls = 0;[mi,n] = size(Aineq);%变量和线性不等式约束的数量F = 0 (1,n);me = size(Aeq,1);%线性等式约束的个数Opts = optimoptions(“linprog”，“算法”，“对偶单纯形”，“显示”，“没有”）;inq_iis = true(mi,1);从问题中的所有不等式开始。Eq_iis = true(me,1);从问题中的所有等式开始。为I =1:mi ineq_iis(I) = 0;去除不等式i[~, ~, exitflag] = linprog (f, Aineq (ineq_iis:), bineq (ineq_iis),.．.Aeq,说真的,磅,乌兰巴托,[],选择);Ncalls = Ncalls + 1;如果Exitflag == 1%如果现在可行Ineq_iis (i) = 1;返回i到问题结束结束为I =1:me eq_iis(I) = 0;%删除等号i[~，~，exitflag] = linprog(f,Aineq,bineq，.．.Aeq (eq_iis:),说真的(eq_iis)磅,乌兰巴托,[],选择);Ncalls = Ncalls + 1;如果Exitflag == 1%如果现在可行Eq_iis (i) = 1;返回i到问题结束结束结束

此代码创建elasticfilterhelper函数。

函数[ineq_iis,eq_iis,ncalls,fval0] = elasticfilter(Aineq,bineq,Aeq,beq,lb,ub) ncalls = 0;[mi,n] = size(Aineq);%变量和线性不等式约束的数量me = size(Aeq,1);Aineq_r = [Aineq -1.0*eye(mi) zero (mi,2*me)];Aeq_r = [Aeq zero (me,mi) eye(me) -1.0*eye(me)];每个等式约束有两条裤子。Lb_r = [lb(:);0 (mi + 2 *我,1)];Ub_r = [ub(:);正(mi + 2 *我,1)];Ineq_slack_offset = n;Eq_pos_slack_offset = n + mi;Eq_neg_slack_offset = n + mi + me;F = [0 (1,n) 1 (1,mi+2*me)];Opts = optimoptions(“linprog”，“算法”，“对偶单纯形”，“显示”，“没有”）;Tol = 1e-10;inq_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;而Exitflag == 1 && fval > tol%可行，一些宽松是非零的C = 0;为I = 1:mi j = ineq_slack_offset+ I;如果X (j) > tol ub_r(j) = 0.0;Ineq_iis (i) = true;C = C +1;结束结束为I = 1:me j = eq_pos_slack_offset+ I;如果X (j) > tol ub_r(j) = 0.0;Eq_iis (i) = true;C = C +1;结束结束为I = 1:me j = eq_neg_slack_offset+ I;如果X (j) > tol ub_r(j) = 0.0;Eq_iis (i) = true;C = C +1;结束结束[x, fval exitflag] = linprog (f Aineq_r bineq Aeq_r,说真的,lb_r, ub_r,[],选择);如果Fval > 0 fval0 = Fval;结束Ncalls = Ncalls + 1;结束结束

此代码创建generatecoverhelper函数。代码使用相同的索引技术跟踪约束删除循环中的IIS代码。

函数[coversetineq,coverseteq,nlp] = generatecover(Aineq,bineq,Aeq,beq,lb,ub)返回线性不等式的封面集，线性%等式，调用linprog的总次数。%改编自Chinneck[1]算法7.3。步数出自这本书。Coversetineq = [];Coverseteq = [];activeA = true(size(bineq));activeAeq = true(size(beq));算法7.3中的步骤1[ineq_iis,eq_iis,ncalls] = elasticfilter(Aineq,bineq,Aeq,beq,lb,ub);NLP = ncalls;Ninf = sum(ineq_iis(:)) + sum(eq_iis(:)));如果Ninf == 1 coversetineq = ineq_iis;Coverseteq = eq_iis;返回结束Holdsetineq = find(ineq_iis);Holdseteq = find(eq_iis);Candidateineq = holdsetineq;Candidateeq = holdseteq;算法7.3中的步骤2而Sum (candidateineq(:)) + Sum (candidateeq(:)) > 0 minsinf = inf;inqflag = 0;为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));如果Fval == 0 coversetineq = [coversetineq;candidateineq(i)];返回结束如果Fval < minsinf ineqflag = 1;Winner = candidateineq(i);Minsinf = fval;Holdsetineq = ineq_iis;如果Numel (ineq_iis(:)) + Numel (eq_iis(:)) == 1 nextwinner = ineq_iis;Nextwinner2 = eq_iis;Nextwinner = [Nextwinner,nextwinner2];其他的Nextwinner = [];结束结束activeA(candidateineq(i)) = true;结束为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));如果Fval == 0 coverseteq = [coverseteq;candidateeq(i)];返回结束如果Fval < minsinf ineqflag = -1;Winner = candidateeq(i);Minsinf = fval;Holdseteq = eq_iis;如果Numel (ineq_iis(:)) + Numel (eq_iis(:)) == 1 nextwinner = ineq_iis;Nextwinner2 = eq_iis;Nextwinner = [Nextwinner,nextwinner2];其他的Nextwinner = [];结束结束activeAeq(candidateeq(i)) = true;结束算法7.3中的步骤3如果Ineqflag == 1 coversetineq = [coversetineq;winner];activeA(winner) = false;如果Nextwinner coversetineq = [coversetineq; Nextwinner];返回结束结束如果Ineqflag == -1 coverseteq = [coverseteq;winner];activeAeq(winner) = false;如果Nextwinner coverseteq = [coverseteq; Nextwinner];返回结束结束Candidateineq = holdsetineq;Candidateeq = holdseteq;结束结束

另请参阅

linprog