ドキュメンテーション センター

  • 評価版
  • 製品アップデート

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

0-1 整数計画

この例では、関数 bintprog を使用した 0-1 整数計画を使って割り当ての問題を解く方法を説明します。

オフィス割り当て問題

Marcelo、Rakesh、Peter、Tom、Marjorie、および Mary Ann の 6 人を 7 つのオフィスに割り当てたいと思います。各オフィスに配置できるのは 1 人だけで、各人は正確に 1 つのオフィスにのみ割り当てられます。これらの人は、オフィスを選ぶことができます。彼らの選好を考慮する際は、年功 (MathWorks に勤務した期間が長いほど、年功が高くなります) が基準とされます。窓のあるオフィスと窓のないオフィスがあります。また、一部の窓は他の窓より大きいです。さらに、Peter と Tom はよく一緒に働いているため、1 つだけ離れたオフィスに配置したいと思います。Marcelo と Rakesh もよく一緒に働いているため、同様にします。

オフィスのレイアウト

オフィス 1、2、3、および 4 は内部オフィス (窓のない) です。オフィス 5、6、および 7 には窓がありますが、オフィス 5 の窓は他の 2 つのオフィスよりも小さいです。ここで、これらのオフィスの配置方法を示します。

text(0.1, .73, 'office1');
text(.35, .73, 'office2');
text(.60, .73, 'office3');
text(.82, .73, 'office4');

text(.35, .42, 'office5');
text(.60, .42, 'office6');
text(.82, .42, 'office7');
title('Office layout: window offices are in the bottom row');
axis off
set(gcf,'color','w');

問題の定式化

問題の定式化を選択しなければなりません。最初の手順では、問題において解変数 x の各要素が表すものを選択します。これは 0-1 整数問題であるため、オフィスに割り当てられる人を各要素が表すようにすることが適切な選択です。社員がオフィスに割り当てられると変数の値が 1 になり、割り当てられないと 0 になります。各人についての検討は、常に以下の順序で行います。

1. Mary Ann
2. Marjorie
3. Tom
4. Peter
5. Marcelo
6. Rakesh

x はベクトルでなければなりません。次に、x(1)x(7) は、オフィス 1 ~オフィス 7 に割り当てられている Mary Ann に相当します。次の 7 つの要素は、オフィス 1 ~オフィス 7 に割り当てられている Marjorie に相当します (以下、他の人についても同様)。6 人を 7 つのオフィスに割り当てるため、x ベクトルの要素は合計で 42 個になります。

年功

MathWorks で勤務している期間が長いほど、選好が重視されるように、年功に基づいて選好に重みを付けたいと思います。年功は以下のようになっています。Mary Ann は 9 年、Marjorie は 10 年、Tom は 5 年、Peter は 3 年、Marcelo は 1.5 年、Rakesh は 2 年です。年功に基づいて、正規化された重みベクトルを作成できます。

seniority = [9 10 5 3 1.5 2];
weightvector = seniority/sum(seniority);

人々のオフィスの好み

行がオフィスに対応し、列が人に対応するような選好行列を設定します。各社員の希望は数値で示されています。第一希望のオフィスの数値が最も大きく、すべての数値の合計、つまりその社員の列の合計が 100 になるようにしています。各人の選好が列ベクトルにリストされます。

MaryAnn = [0; 0; 0; 0; 10; 40; 50];
Marjorie = [0; 0; 0; 0; 20; 40; 40];
Tom = [0; 0; 0; 0; 30; 40; 30];
Peter = [1; 3; 3; 3; 10; 40; 40];
Marcelo = [3; 4; 1; 2; 10; 40; 40];
Rakesh = [10; 10; 10; 10; 20; 20; 20];

各人の選好ベクトルの i 番目の要素は、その人物が i 番目のオフィスをどの程度評価しているかを示しています。つまり、選好行列を組み合わせたものは次のようになります。

prefmatrix = [MaryAnn Marjorie Tom Peter Marcelo Rakesh];

年功が列の基準となるように、weightvector を使用して選好行列に重みを付けたいと思います。また、この行列を列順にベクトルに変換して、x ベクトルに対応するようにすれば、より都合が良いでしょう。

PM = prefmatrix * diag(weightvector);
c = PM(:);

目的関数

目的は、年功によって重みが付けられた選好の満足度を最大限にすることです。以下は線形目的関数です。

     max c'*x

これは、次の関数と等価です。

     min -c'*x

制約

最初の制約セットでは、各人が正確に 1 つのオフィスを得ることが求められます。つまり、各人について、該当する人に対応する x 値の合計が正確に 1 であることが求められます。たとえば、2 番目の人である Marjorie の場合、これは sum(x(8:14))=1 であることを意味します。適切な行列を作成することで、これらの線形制約を等式行列 Aeq と右側のベクトル beq (ここで、Aeq*x = beq) で表すことができます。行列 Aeq は、1 と 0 で構成されます。たとえば、Aeq の 2 番目の行は、正確に 1 つのオフィスを得る Marjorie に対応しており、行は以下のようになります。

     0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 ... 0 0 0

つまり、列 8 ~ 14 に 7 個の 1 があり、他の列は ゼロです。Aeq(2,:)*x = 1 は、sum(x(8:14)) = 1 と等価です。

numOffices = 7;
numPeople = 6;
numDim = numOffices * numPeople;
onesvector = ones(1,numOffices);
% Each row of Aeq corresponds to one person.
Aeq = blkdiag(onesvector,onesvector,onesvector,onesvector, ...
    onesvector,onesvector);
beq = ones(numPeople,1);
% View the structure of Aeq, that is, where there are nonzeros (ones)
figure;
spy(Aeq)
set(gcf,'color','w');
set(get(gcf,'CurrentAxes'),...
    'PlotBoxAspectRatio',[43 8 1],...
    'YTick',[1 2 3 4 5 6],...
    'YTickLabel',{'Mary Ann','Marjorie','Tom','Peter','Marcelo','Rakesh'});
title('Equality constraints: each person gets exactly one office')

第 2 の制約セットは不等式です。これらの制約では、各オフィスに 1 人だけ配置される (つまり、各オフィスに 1 人がいるか誰もいない) ように指定されます。これらの制約を表現するために、A*x <= b となるような行列 A とベクトル b を作成します。行列 Ab の各行は、オフィスに対応します。したがって、行 1 はオフィス 1 に割り当てられる社員を示します。どの行も一定のパターンで 1 と 0 が並んでいます。行 1 は次のようになります。

    1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0  ... 1 0 0 0 0 0 0

この後の各行は、これと似ていますが、右に 1 つシフト (循環的に) したものになります。たとえば、行 3 はオフィス 3 に対応しており、A(3,:)*x <= 1 (つまり、オフィス 3 には 1 人しか配置できないこと) を示しています。

A = repmat(eye(numOffices),1,numPeople);
b = ones(numOffices,1);
% View the structure of Aeq, that is, where there are nonzeros (ones)
figure;
spy(A)
set(gcf,'color','w');
set(get(gcf,'CurrentAxes'),...
    'YTick',[1 3 5 7],...
    'YTickLabel',{'Office 1','Office 3','Office 5','Office 7'},...
    'XTick',[1 8 15 22 29 36],...
    'XTickLabel',{'Mary Ann','Marjorie','Tom','Peter','Marcelo','Rakesh'});
title('Inequality constraints: no more than one person per office')

次の制約セットも不等式であるため、それらを、上記から既に不等式を含んでいる行列 A とベクトル b に追加します。Tom と Peter には、1 つだけ離れたオフィスを割り当てたいと思います。Marcelo と Rakesh についても同様です。まず、彼らの物理的位置に基づき、近似マンハッタン距離を使用して、オフィスの距離行列を作成します。これは対称行列です。

D = zeros(numOffices);
%   Set up the top right half of the matrix
D(1,2:end) = [1 2 3 2 3 4];
D(2,3:end) = [1 2 1 2 3];
D(3,4:end) = [1 2 1 2];
D(4,5:end) = [3 2 1];
D(5,6:end) = [1 2];
D(6,end)   = 1;
% The lower left half is the same as the upper right
D = triu(D)' + D;

% We find the offices that are more than one distance unit away.
[officeA,officeB] = find(D > 1);
numPairs = length(officeA)
numPairs =

    26

これにより、オフィスの numPairs ペアが見つかります。これらの numPairs に関して、Tom がペアの一方のオフィスを占有している場合、Peter はペアの他方のオフィスを占有できません。そうでなければ、彼らは、2 つ以上離れたオフィスに配置されることになってしまいます。Marcelo と Rakesh についてもこれと同様です。これにより、さらに 2*numPairs 多い不等式制約が与えられ、これらの制約を Ab に追加することになります。

% Add enough rows to A to accommodate these constraints.
numrows = 2*numPairs + numOffices;
A((numOffices+1):numrows, 1:numDim) = zeros(2*numPairs,numDim);

numPairs の各オフィス ペア、officeA の Tom に対応する x(i)、および officeB の Peter に対応する x(j) について、制約は以下のようになります。

x(i) + x(j) <= 1

つまり、Tom または Peter のいずれかがこれらのオフィスの 1 つを占有できますが、両人とも占有することはできません。

% Tom is person 3 and Peter is person 4
tom = 3;
peter = 4;
% Marcelo is person 5 and Rakesh is person 6
marcelo = 5;
rakesh = 6;
% The following anonymous functions return the index in x
% corresponding to Tom, Peter, Marcelo and Rakesh respectively
% in office i.
tom_index=@(officenum) (tom-1)*numOffices+officenum;
peter_index=@(officenum) (peter-1)*numOffices+officenum;
marcelo_index=@(officenum) (marcelo-1)*numOffices+officenum;
rakesh_index=@(officenum) (rakesh-1)*numOffices+officenum;

for i = 1:numPairs
    tomInOfficeA = tom_index(officeA(i));
    peterInOfficeB = peter_index(officeB(i));
    A(i+numOffices, [tomInOfficeA, peterInOfficeB]) = 1;

    % Repeat for Marcelo and Rakesh, adding more rows to A and b.
    marceloInOfficeA = marcelo_index(officeA(i));
    rakeshInOfficeB = rakesh_index(officeB(i));
    A(i+numPairs+numOffices, [marceloInOfficeA, rakeshInOfficeB]) = 1;
end
b(numOffices+1:numOffices+2*numPairs) = ones(2*numPairs,1);
% View the structure of the newly added constraints in A, that is,
%   where there are nonzeros (ones)
figure;
spy( A((numOffices+1):numrows,:) )
set(gcf,'color','w');
title('Inequality constraints: Tom & Peter nearby; Marcelo & Rakesh nearby')

BINTPROG を使用した解法

問題の定式化は、以下の線形目的関数で構成されます。

   min -c'*x

これは、以下の線形制約に従います。

   Aeq*x = beq
   A*x <= b

これは BINTPROG で期待される形式であるため、これらの行列を BINTPROG に渡すことができます。

% Let BINTPROG choose the start point.
x0 = [];
f = -c;
% Show the iterative output for each node displayed in the branch and
% bound algorithm.
options = optimoptions('bintprog','Display','iter','NodeDisplayInterval',1);
[x,fval,exitflag,output] = bintprog(f,A,b,Aeq,beq,x0,options);
fval
exitflag
output
% View the solution to see who got what office.
printofficeassign(x);
title('Solution for default BranchStrategy and NodeSearchStrategy');
Explored   Obj of LP   Obj of best   Unexplored   Best lower    Relative gap
 nodes    relaxation  integer point     nodes    bound on obj  between bounds
      1       -33.87            -          2          -33.87          -
      2       -32.48            -          3          -33.87          -
*     3       -33.84       -33.84          0          -33.87       0.094%
Optimization terminated.

fval =

  -33.8361


exitflag =

     1


output = 

          iterations: 20
               nodes: 3
                time: 0.6552
           algorithm: 'LP-based branch-and-bound'
      branchStrategy: 'maximum integer infeasibility'
    nodeSrchStrategy: 'best node search'
             message: 'Optimization terminated.'

この問題の場合、年功による選好の満足度が -fval の値まで最大化されます。正の exitflag は、BINTPROG が収束したことを示します。この出力構造体により、調べられたノードの数、計算にかかった時間、および LP 緩和サブ問題を解く際に使用された累積反復の回数に関する情報を得ることができます。

オプションの変更

調べられるノードの数、計算にかかる時間、または反復の回数を低減するために、1 組のオプションの変更を試みることができます。BINTPROG では分岐限定アルゴリズムが使用されていますが、オプションを使用して、このアルゴリズムを調整できます。たとえば、既定の分岐法 [maxinfeas] では、次に分岐することが不可能な最大の整数の変数 (つまり、値が 0.5 に最も近い変数) を選択します。この分岐法を [mininfeas] に設定して問題を再び実行できます。[mininfeas] では、次に分岐することが不可能な最小の整数の変数 (つまり、値が 0 または 1 に最も近いが 0 または 1 に等しくない変数) を選択します。

% Try BranchStrategy = mininfeas
options.BranchStrategy = 'mininfeas';
[x,fval,exitflag,output] = bintprog(f,A,b,Aeq,beq,x0,options);
fval
exitflag
output
printofficeassign(x);
title('Solution for BranchStrategy=mininfeas and default NodeSearchStrategy');
Explored   Obj of LP   Obj of best   Unexplored   Best lower    Relative gap
 nodes    relaxation  integer point     nodes    bound on obj  between bounds
      1       -33.87            -          2          -33.87          -
*     2       -33.79       -33.79          1          -33.87        0.24%
      3       -33.86       -33.79          2          -33.87        0.24%
*     4       -33.84       -33.84          1          -33.86       0.071%
Optimization terminated.

fval =

  -33.8361


exitflag =

     1


output = 

          iterations: 21
               nodes: 4
                time: 0.1092
           algorithm: 'LP-based branch-and-bound'
      branchStrategy: 'minimum integer infeasibility'
    nodeSrchStrategy: 'best node search'
             message: 'Optimization terminated.'

この問題の場合、代替の分岐法によってノード数と反復回数が減らされますが、以前と同じ解は見つけられます。

最後に、別のノード探索法を選択できます。これは、探索木で次の探索対象ノードを選択するために分岐限定アルゴリズムで使用されるノード探索法です。既定の設定では、最良ノード探索法 "bn" を使用して探索が行われます。これは、次に探索される、目的関数の最小の下限をもつノードを選択する探索法です。深度探索法 "df" を使用するように、これを変更できます。探索木の各ノードで、木の 1 レベル下にまだ探索されていない子ノードが存在する場合、このアルゴリズムは、その子を選択して探索します。そうでない場合、このアルゴリズムは、木の 1 レベル上のノードに移動し、そのノードから 1 レベル下の子ノードを選択します。

% Try NodeSearchStrategy = df
options.NodeSearchStrategy = 'df';
[x,fval,exitflag,output] = bintprog(f,A,b,Aeq,beq,x0,options);
fval
exitflag
output
printofficeassign(x);
title('Solution for BranchStrategy=mininfeas and NodeSearchStrategy=df');
Explored   Obj of LP   Obj of best   Unexplored   Best lower    Relative gap
 nodes    relaxation  integer point     nodes    bound on obj  between bounds
      1       -33.87            -          2          -33.87          -
      2       -33.86            -          3          -33.87          -
      3       -33.84            -          4          -33.87          -
      4       -33.24            -          5          -33.87          -
      5       -32.99            -          6          -33.87          -
      6       -32.99            -          7          -33.87          -
      7       -32.79            -          8          -33.87          -
      8       -32.41            -          9          -33.87          -
      9       -30.29            -         10          -33.87          -
     10       -30.29            -         11          -33.87          -
     11       -29.74            -         12          -33.87          -
     12       -29.74            -         13          -33.87          -
*    13       -29.02       -29.02         12          -33.87          16%
     14       -22.82       -29.02         11          -33.87          16%
     15       -22.82       -29.02         10          -33.87          16%
     16       -24.46       -29.02          9          -33.87          16%
     17       -24.46       -29.02          8          -33.87          16%
     18       -29.22       -29.02          9          -33.87          16%
     19       -22.99       -29.02          8          -33.87          16%
     20       -24.03       -29.02          7          -33.87          16%
     21        -26.1       -29.02          6          -33.87          16%
*    22       -29.67       -29.67          5          -33.87          14%
     23       -29.42       -29.67          4          -33.87          14%
     24       -29.67       -29.67          3          -33.87          14%
     25       -32.48       -29.67          4          -33.87          14%
     26       -31.77       -29.67          5          -33.87          14%
     27       -31.77       -29.67          6          -33.87          14%
     28       -28.97       -29.67          5          -33.87          14%
     29       -26.05       -29.67          4          -33.87          14%
     30       -28.97       -29.67          3          -33.87          14%
     31       -28.97       -29.67          2          -33.87          14%
*    32       -33.84       -33.84          1          -33.87       0.094%
Optimization terminated.

fval =

  -33.8361


exitflag =

     1


output = 

          iterations: 141
               nodes: 32
                time: 0.2964
           algorithm: 'LP-based branch-and-bound'
      branchStrategy: 'minimum integer infeasibility'
    nodeSrchStrategy: 'depth first search'
             message: 'Optimization terminated.'

この問題の場合、代替のノード探索法を使用して、ノード数、反復回数、および時間が低減されます。以前と同じ解が見つけられます。

この情報は役に立ちましたか?