Main Content

線形実行不可能性の調査

ライブ スクリプトを開く

この例では、問題が実行不可能となる線形制約を調査する方法を説明します。これらの手法の詳細については、Chinneck [1] および [2] を参照してください。

線形制約によって問題が実行不可能になる場合、サブセットの一部のメンバーを削除すれば残りのメンバーは問題が実行可能になるという条件で、実行不可能な制約のサブセットを探したい場合があります。そのようなサブセットを、"制約の既約実行不可能サブセット" (IIS) といいます。反対に、実行可能な制約の最大濃度サブセットを見つける必要がある場合もあります。このサブセットは、"最大実行可能サブセット" (MaxFS) と呼ばれます。この 2 つの概念は関連していますが、同一ではありません。問題によっては、多数の異なる IIS が含まれ、濃度も異なる場合があります。

この例では、IIS を見つける 2 つの方法と、実行可能な制約のセットを取得する 2 つの方法を説明します。この例では、指定されたすべての境界が正しい、つまり、lb 引数と ub 引数にエラーはないことを前提とします。

実行不可能な例

サイズが 150 行 15 列の線形不等式を表す乱数行列 A を作成します。対応するベクトル 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 から N までの各 i について次を行います。

  • 問題から制約 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

境界制約とともに問題を実行不可能にしている不等式制約は 1 つだけです。制約は次のとおりです。

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

境界とともにこの制約が実行不可能なのはなぜでしょうか。A の当該行について、その絶対値の和を求めます。

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

境界があるので、x ベクトルは -1 と 1 の間の値をとるため、A(iis,:)*x は b(iis) = -10 未満にはなり得ません。

deletionfilter は linprog 呼び出しを何回実行したでしょうか。

disp(ncall)
   150

問題には 150 件の線形制約があるため、この関数は linprog を 150 回呼び出しました。

弾性フィルター

すべての制約を検証する削除フィルターの代替として、弾性フィルターを試してみましょう。このフィルターは、次のように機能します。

まず、各制約 i が正の量 e(i) によって違反されることを許容します。ここで、等式制約は、加法的な正の弾性値と減法的な正の弾性値をもちます。

Aineqx≤bineq+eAeqx=beq+e1-e2

次に、以下の関連する線形計画問題 (LP) を解きます。

minx,e∑ei

これは、リストされた制約および ei≥0 に従います。

  • 関連する LP に解があれば、厳密に正である関連 ei をもつすべての制約を削除し、それらの制約をインデックスのリストに記録します (潜在的な IIS 要素)。前の手順に戻り、簡略化された新しい関連する LP を解きます。

  • 関連する LP に解がない (実行不可能である)、または厳密に正である関連 ei がない場合は、フィルターを終了します。

弾性フィルターは、削除フィルターよりも大幅に少ない反復回数で終了できます。弾性フィルターは一度に多くのインデックスを 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 では 5 つの制約がリストされますが、削除フィルターでは 1 つの制約のみ検出されます。返されたセットに対して削除フィルターを実行して、本当の 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

弾性フィルターの結果における 5 番目の制約、不等式 114 がその IIS です。この結果は、削除フィルターの回答とも一致します。これらの手法の違いは、弾性フィルターと削除フィルターを組み合わせた場合に、使用される linprog 呼び出しの回数が大幅に少なくなるということです。弾性フィルターによって使用され、次いで削除フィルターによって使用される linprog 呼び出しの合計回数を表示します。

disp(ncall + ncall2)
     7

ループでの IIS の削除

一般に、1 つの IIS を取得するだけでは、最適化問題が解けない理由のすべてを見つけることはできません。実行不可能な問題を修正するには、IIS を見つけて問題から削除することを、その問題が実行可能になるまで繰り返します。

次のコードは、問題が実行可能になるまで、その問題から一度に 1 つずつ 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 を同時に削除していることから、これら 2 つの制約が IIS となることがわかります。

MaxFS の検出

実行可能な制約のセットを取得するためのもう 1 つの手法は、MaxFS を直接見つけることです。Chinneck [1] で説明されているように、MaxFS の検出は NP 完全問題であり、問題が MaxFS を見つけるための効率的なアルゴリズムを必ずしももつわけではないことを意味します。しかし、Chinneck は、効率的に機能し得るいくつかのアルゴリズムを推奨しています。

Chinneck の Algorithm 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.

このカバー セットが弾性フィルターで得られた iiselastic セットと完全に同一であることに注意してください。一般に、弾性フィルターが検出するカバー セットは大きすぎます。Chinneck の Algorithm 7.3 は、弾性フィルターの結果を使用して開始し、必要な制約のみを保持します。

Chinneck の Algorithm 7.3 は、linprog に対する呼び出しを 40 回実行して MaxFS の計算を完了します。この回数は、ループで IIS を削除するプロセスで先に使用した 28 回という呼び出し回数よりもいくらか多くなっています。

ループで削除された不等式は、Algorithm 7.3 によって削除された不等式と完全に同一ではないことにも注意してください。ループでは不等式 114、97、82、60、および 64 を削除する一方、Algorithm 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.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

参考

関連するトピック