调查线性Infeasibilities

打开实时脚本

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

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

此示例示出了找到一个IIS的两种方式，并获得一个可行的一组约束的方法有两种。这个例子假定所有给定的边界是正确的，这意味着磅和乌兰巴托论据有没有错误。

不可行例

创建一个随机矩阵一种表示大小为150 × 15的线性不等式。设置相应的向量B.将这些值的5%更改为–10。

N=15；rng违约A = randn（[10 * N，N]）;B = 10吨*酮（尺寸（A，1），1）;AEQ = [];BEQ = [];B（兰特（大小（B））<= 0.05）= -10;F =酮（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辅助函数在这个例子到此结束。

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

见效果deletionfilter上示例性数据。

[ineqs、eqs、ncall]=删除筛选器（A、b、Aeq、beq、lb、ub）；

这个问题没有等式约束。找到不等式约束的指数和b（iis）。

IIS =查找（ineqs）

IIS = 114

b（iis）

ans = -10

只有一个不等式约束和有界约束导致问题不可行。约束条件是

A（iis，：）*x<=b（iis）。

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

DISP（总和（ABS（A（IIS，:)）））

8.4864

由于边界的限制X矢量-1之间具有值和1，依此A（iis，：）*x不能少于b（iis）= –10.

多少linprog电话打了deletionfilter表演

显示（ncall）

这个问题有150个线性约束，所以称为功能linprog150次。

弹性过滤器

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

首先，允许每个约束一世被正数违反的E（i），其中等式约束同时具有加法和减法正弹性值。

$\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®代码实现此过滤器，请参阅弹性过滤器辅助函数在这个例子到此结束。

见效果弹性过滤器上示例性数据。

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

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

iiselastic=查找（ineqselastic）

iiselastic =5×12 60 82 97 114

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

A_1=A（ineqseastic>0，：）；b_1=b（ineqseastic>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弹性筛选器使用的调用，后跟删除筛选器。

显示（ncall+ncall 2）

7.

在循环中删除IIS

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

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

代码使用布尔向量跟踪问题中的原始变量活跃的表示当前约束（行）的一种矩阵和布尔向量activeAeq来表示当前的约束AEQ矩阵。添加或移除约束时，代码索引为一种或AEQ从而使行数不改变，即使约束的数量变化。

运行此代码将返回idx2，中非零元素指数的向量活跃的：

idx2 =找到(activeA)

假设变量布尔向量的长度和idx2. 然后

IDX2（FIND（VAR））

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

OPTS = optimoptions（'linprog'那“显示”那“没有”); activeA=真（大小（b））；activeAeq=真（尺寸（beq））；[~，~，exitflag]=linprog（f，A，b，Aeq，beq，lb，ub，opts）；ncl=1；虽然Exitflag < 0 [ineqselastic,eqselastic,ncall] =...elasticfilter（A（activeA，:)，B（activeA），AEQ（activeAeq，:)，BEQ（activeAeq），LB，UB）;NCL = NCL + ncall;idxaa =查找（activeA）;idxae =找到（activeAeq）;TMPA = idxaa（FIND（ineqselastic））;tmpae = idxae（找到（eqselastic））;AA = A（TMPA，:);BB = B（TMPA）;AE = AEQ（tmpae，:);是= BEQ（tmpae）; [ineqs,eqs,ncall] =...deletionfilter (AA、bb AE,磅,乌兰巴托);NCL = NCL + ncall;activeA (tmpa (ineqs)) = false;activeAeq (tmpae(方程式))= false;disp ([“删除不平等”，int2str（（TMPA（ineqs））'），“和平等”，int2str（（tmpae（eqs）））][~，~，exitflag]=...linprog（f，A（activeA，：），b（activeA），Aeq（activeAeq，：），beq（activeAeq），lb，ub，opts）；ncl=ncl+1；结尾

删除不等式114和等式删除不等式97和等式删除不等式64 82和等式删除不等式60和等式

fprintf（'linprog调用数：%d\n'，ncl）

linprog呼叫数：28

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

找到maxf

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

使用Chinneck's算法7.3求a封面套装当移除时，给出可行集的约束。该算法在计算机上实现发电机换向辅助函数在这个例子到此结束。

[coversetineq，coverseteq，NLP] = generatecover（A，B，AEQ，BEQ，LB，UB）

coversetineq =5×1114 97 60 82 2

coverseteq = []

nlp=40

删除这些约束和解决LP。

usemeineq =真（大小（B））;usemeineq（coversetineq）= FALSE;去掉不等式约束usemeeq=真（尺寸（beq））；usemeeq（coverseteq）=假；％删除等式约束[x, fvals exitflags] =...linprog（F，A（usemeineq，:)，B（usemeineq），AEQ（usemeeq），BEQ（usemeeq），LB，UB）;

找到了最优解。

请注意，封面集与封面集完全相同自言自语的从设置弹性过滤器。通常，弹性过滤器会发现太大的覆盖集。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(大小(b));usemeineq([2] 82年)= true;ineqs = deletionfilter (A (usemeineq:)、b (usemeineq) Aeq,说真的,磅,乌兰巴托);disp (ineqs)

1 1

参考文献

[1] 金内克，J.W。可行性和不可行性的优化：算法和计算方法。斯普林格，2008年。

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

辅助函数

此代码创建deletionfilter助手功能。

作用[ineq_iis，eq_iis，ncalls] = deletionfilter（Aineq，bineq，AEQ，BEQ，LB，UB）ncalls = 0;[MI中，n] =尺寸（Aineq）;%变量数与线性不等式约束f=零（1，n）；me=尺寸（Aeq，1）；%线性等式约束数OPTS = optimoptions（“linprog”那“算法”那“双单”那“展示”那“没有”); ineq_iis=真（mi，1）；%从问题中的所有不等式开始eq_iis=真（me，1）；%从问题中的所有等式开始为了i=1:mi ineq_iis（i）=0；％删除我的不平等[~，~，exitflag]=linprog（f，Aineq（ineq_-iis，：），bineq（ineq_-iis），...Aeq、beq、lb、ub、[]和opts）；ncalls=ncalls+1；如果exitflag = = 1％如现在可行ineq_iis（1）= 1;返回i到问题结尾结尾为了I =1: eq_iis(I) = 0;删除等号i[〜，〜，exitflag] = linprog（F，Aineq，bineq，...AEQ（eq_iis，:)，BEQ（eq_iis），LB，UB，[]，OPTS）;ncalls = ncalls + 1;如果exitflag = = 1％如现在可行等式（i）=1；返回i到问题结尾结尾结尾

此代码创建弹性过滤器助手功能。

作用[ineq_iis，eq_iis，ncalls，fval0] = elasticfilter（Aineq，bineq，AEQ，BEQ，LB，UB）ncalls = 0;[MI中，n] =尺寸（Aineq）;%变量数与线性不等式约束me=尺寸（Aeq，1）；Aineq_r=[Aineq-1.0*眼睛（mi，2*me）]；Aeq_r=[Aeq零（me，mi）眼睛（me）-1.0*眼睛（me）]；%每个相等约束有两条松弛lb_r=[lb（：）；零（mi+2*me，1）]；ub_r=[ub（：）；inf（mi+2*me，1）]；ineq_松弛_偏移=n；等式位置松弛偏移=n+mi；等式负松弛偏移=n+mi+me；f=[零（1，n）一（1，mi+2*me）]；opts=options(“linprog”那“算法”那“双单”那“展示”那“没有”）;托尔= 1平台以及;ineq_iis = false (mi, 1);eq_iis = false(我,1);[x, fval exitflag] = linprog (f Aineq_r bineq Aeq_r,说真的,lb_r, ub_r,[],选择);fval0 = fval;ncalls = ncalls + 1;虽然exitflag==1&&fval>tol％可行和一些松弛是非零C = 0;为了i=1:mi j=ineq\u slack\u offset+i；如果x（j）>tol ub_r（j）=0.0；ineq_iis（i）=真；c=c+1；结尾结尾为了i=1:me j=eq\u pos\u slack\u offset+i；如果x（j）>tol ub_r（j）=0.0；等式（i）=真；c=c+1；结尾结尾为了i=1:me j=eq\u neg\u slack\u offset+i；如果x（j）>tol ub_r（j）=0.0；等式（i）=真；c=c+1；结尾结尾[x，fval，exitflag]=linprog（f，Aineq_r，bineq，Aeq_r，beq，lb_r，ub_r，[]，opts）；如果fval>0 fval0=fval；结尾ncalls = ncalls + 1;结尾结尾

此代码创建发电机换向辅助函数。该代码使用了与代码相同的索引技术来跟踪约束在循环中删除IIS密码

作用[coversetineq，coverseteq，NLP] = generatecover（Aineq，bineq，AEQ，BEQ，LB，UB）%返回线性不等式的覆盖集，线性不等式的覆盖集％等式，并呼吁linprog的总数。％从Chinneck [1] 7.3算法改编。步骤编号从本书。coversetineq=[]；coverseteq=[]；activeA=true（大小（bineq））；activeAeq=true（大小（beq））；%算法7.3的步骤1[ineq_-iis，eq_-iis，ncalls]=弹性过滤器（Aineq，bineq，Aeq，beq，lb，ub）；nlp=ncalls；ninf=sum（ineq_-iis（：）+sum（eq_-iis（：）；如果Ninf == 1 coversetineq = inq_iis;coverseteq = eq_iis;回来结尾holdsetineq=find（ineq_-iis）；holdseteq=find（eq_-iis）；candidateineq=holdsetineq；候选资格Q=holdseteq；％7.3算法的步骤2中虽然sum（candidateq（：）+sum（candidateq（：）>0分钟inf=inf；ineqflag=0；为了i = 1:length(candidateineq(:)) activeA(candidateineq(i)) = false;idx2 =找到(activeA);idx2eq =找到(activeAeq);[ineq_iis, eq_iis ncalls fval] = elasticfilter (Aineq (activeA:), bineq (activeA) Aeq (activeAeq),说真的(activeAeq)磅,乌兰巴托);NLP = NLP + ncalls;ineq_iis = idx2(找到(ineq_iis));eq_iis = idx2eq(找到(eq_iis));如果FVAL == 0 coversetineq = [coversetineq; candidateineq（I）];回来结尾如果Fval < minsinf inqflag = 1;赢家= candidateineq(我);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（candidateq（：））activeAeq（candidateq（i））=false；idx2=查找（activeA）；idx2eq=查找（activeAeq）；[ineq_iis，eq_iis，NCALL，fval]=弹性过滤器（Aineq（activeA），bineq（activeA），Aeq（activeAeq，：），beq（activeAeq），lb，ub）；nlp=nlp+NCALL；ineq_iis=idx2（查找（ineq_iis））；eq_iis=idx2eq（find（eq_iis））；如果fval==0 coverseteq=[coverseteq；候选资格（i）]；回来结尾如果fval如果numel（ineq_iis（：）+numel（eq_iis（：）==1 nextwinner=ineq_iis；nextwinner2=eq_iis；nextwinner=[nextwinner，nextwinner2]；其他的nextwinner = [];结尾结尾activeAeq（候选资格（i））=真；结尾%算法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;结尾结尾