调查线性不可行性

打开实时脚本

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

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

这个示例展示了查找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路径上。本例后面使用的其他辅助函数也是如此。

看看效果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

由于边界的关系x向量的值在-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(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控件的当前约束(行)一个矩阵和布尔向量activeAeq的当前约束Aeq矩阵。当添加或删除约束时，代码索引到一个或Aeq这样行号就不会改变，即使约束的数量改变了。

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

idx2 = find(activeA)

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

idx2(找到(var))

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

Opts = optimoptions(“linprog”，“显示”，“没有”）;activeA = true(大小(b));activeAeq = true(大小(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));Tmpae = idxae(find(eqselastic));AA = A(tmpa，:);Bb = b(tmpa);AE = Aeq(tmpae，:);Be = beq(tmpae); [ineqs,eqs,ncall] =.．.deletionfilter (AA、bb AE,磅,乌兰巴托);NCL = NCL + ncall;activeA(tmpa(ineqs)) = false;activeAeq(tmpae(eqs)) = false;disp ([“消除不平等”, int2str ((tmpa (ineqs))”),“和等式”，int2str((tmpae(eqs))')]) [~，~，exitflag] =.．.linprog (f (activeA:)、b (activeA) Aeq (activeAeq:),说真的(activeAeq)磅,乌兰巴托,选择);NCL = NCL + 1;结束

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

流(linprog调用数:%d\nncl)

linprog调用数:28

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

找到maxf

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

使用Chinneck's Algorithm 7.3找到一个覆盖集当去掉约束条件时，给出一个可行集。该算法在generatecover的辅助函数。本例结束．

[coverseineq,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(大小(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)

参考文献

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

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”，“算法”，“对偶单纯形”，“显示”，“没有”）;Ineq_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;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;而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:我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代码。

函数[coverseineq,coverseteq,nlp] = generatcover (Aineq,bineq,Aeq,beq,lb,ub)返回线性不等式的覆盖集，线性不等式的覆盖集% equals，以及调用linprog的总次数。%改编自Chinneck[1]算法7.3。步骤编号来自这本书。Coversetineq = [];Coverseteq = [];activeA = true(大小(bineq));activeAeq = true(大小(beq));7.3算法的步骤1[ineq_iis,eq_iis,ncalls] = elasticfilter(Aineq,bineq,Aeq,beq,lb,ub);NLP = ncalls;nif = sum(ineq_iis(:)) + sum(eq_iis(:));如果nif == 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;Ineqflag = 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(查找(Ineq_iis));Eq_iis = idx2eq(查找(Eq_iis));如果Fval == 0 coversetineq = [coversetineq;candidateineq(i)];返回结束如果Fval < minsinf ineqflag = 1;Winner = candidateineq(i);Minsinf = fval;Holdsetineq = ineq_iis;如果Numel (ineq_iis(:)) + nummel (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(查找(Ineq_iis));Eq_iis = idx2eq(查找(Eq_iis));如果Fval == 0 coverseteq = [coverseteq;candidateeq(i)];返回结束如果Fval < minsinf ineqflag = -1;Winner = candidateeq(i);Minsinf = fval;Holdseteq = eq_iis;如果Numel (ineq_iis(:)) + nummel (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;如果next twinner coversetineq = [coversetineq;next twinner];返回结束结束如果Ineqflag == -1 coverseteq = [coverseteq;winner];activeAeq(winner) = false;如果next twinner coverseteq = [coverseteq;next twinner];返回结束结束Candidateineq = holdsetineq;Candidateeq = holdseteq;结束结束

另请参阅

linprog