Main Content

coneprog

2 次錐計画法ソルバー

R2020b 以降

説明

coneprog 関数は、以下で指定された問題の最小値を見つける 2 次錐計画法ソルバーです。

minxfTx

以下の制約に従います。

Asc(i)xbsc(i)dscT(i)xγ(i)AxbAeqx=beqlbxub.

f、x、b、beq、lb、ub はベクトル、A と Aeq は行列です。i ごとに、行列 Asc(i)、ベクトルの dsc(i) と bsc(i)、およびスカラー γ(i) が、secondordercone を使用して作成された 2 次錐制約に含まれています。

錐制約の詳細については、2 次錐制約を参照してください。

x = coneprog(f,socConstraints) は、次のようにエンコードされた socConstraints 内の制約を使用して 2 次錐計画問題を解きます。

  • Asc(i) = socConstraints(i).A

  • bsc(i) = socConstraints(i).b

  • dsc(i) = socConstraints(i).d

  • γ(i) = socConstraints(i).gamma

x = coneprog(f,socConstraints,A,b,Aeq,beq) は、不等式制約 A*x b と等式制約 Aeq*x = beq が課された問題を解きます。不等式が存在しない場合は A = [] および b = [] と設定してください。

x = coneprog(f,socConstraints,A,b,Aeq,beq,lb,ub) は、解が常に lb ≤ x ≤ ub の範囲に入るように、設計変数 x の上限と下限のセットを定義します。等式が存在しない場合には Aeq = [] および beq = [] と設定してください。

x = coneprog(f,socConstraints,A,b,Aeq,beq,lb,ub,options) は、options で指定された最適化オプションを使用して最小化します。optimoptions を使用してこれらのオプションを設定してください。

x = coneprog(problem) は、problem で説明されている構造体 problem の最小値を求めます。

また、[x,fval] = coneprog(___) は、前の構文の入力引数の組み合わせのいずれかを使用して、解 fval = f'*x における目的関数値を返します。

[x,fval,exitflag,output] = coneprog(___) は上記に加え、終了条件を記述する値 exitflag、および最適化プロセスに関する情報を含む構造体 output を返します。

加えて、[x,fval,exitflag,output,lambda] = coneprog(___) は、解 x でフィールドに双対変数が含まれている構造体 lambda を返します。

すべて折りたたむ

1 つの 2 次錐制約を含む問題を設定するには、2 次錐制約オブジェクトを作成します。

A = diag([1,1/2,0]);
b = zeros(3,1);
d = [0;0;1];
gamma = 0;
socConstraints = secondordercone(A,b,d,gamma);

目的関数ベクトルを作成します。

f = [-1,-2,0];

この問題には、線形制約がありません。これらの制約用の空行列を作成します。

Aineq = [];
bineq = [];
Aeq = [];
beq = [];

x(3) の上限と下限を設定します。

lb = [-Inf,-Inf,0];
ub = [Inf,Inf,2];

coneprog 関数を使用して問題を解きます。

[x,fval] = coneprog(f,socConstraints,Aineq,bineq,Aeq,beq,lb,ub)
Optimal solution found.
x = 3×1

    0.4851
    3.8806
    2.0000

fval = -8.2462

解の要素 x(3) は上限値です。錐制約は次の解においてアクティブです。

norm(A*x-b) - d'*x % Near 0 when the constraint is active
ans = -2.5677e-08

複数の 2 次錐制約を含む問題を設定するには、制約オブジェクトの配列を作成します。時間とメモリを節約するために、最初に、最もインデックスの高い制約を作成します。

A = diag([1,2,0]);
b = zeros(3,1);
d = [0;0;1];
gamma = -1;
socConstraints(3) = secondordercone(A,b,d,gamma);

A = diag([3,0,1]);
d = [0;1;0];
socConstraints(2) = secondordercone(A,b,d,gamma);

A = diag([0;1/2;1/2]);
d = [1;0;0];
socConstraints(1) = secondordercone(A,b,d,gamma);

線形目的関数ベクトルを作成します。

f = [-1;-2;-4];

coneprog 関数を使用して問題を解きます。

[x,fval] = coneprog(f,socConstraints)
Optimal solution found.
x = 3×1

    0.4238
    1.6477
    2.3225

fval = -13.0089

目的関数ベクトルと単一の 2 次錐制約を指定します。

f = [-4;-9;-2];
Asc = diag([1,4,0]);
b = [0;0;0];
d = [0;0;1];
gamma = 0;
socConstraints = secondordercone(Asc,b,d,gamma);

線形不等式制約を指定します。

A = [1/4,1/9,1];
b = 5;

問題を解きます。

[x,fval] = coneprog(f,socConstraints,A,b)
Optimal solution found.
x = 3×1

    3.2304
    0.6398
    4.1213

fval = -26.9225

coneprog ソルバーの反復を観察するには、Display オプションを 'iter' に設定します。

options = optimoptions('coneprog','Display','iter');

2 次錐計画問題を作成し、それを options を使用して解きます。

Asc = diag([1,1/2,0]);
b = zeros(3,1);
d = [0;0;1];
gamma = 0;
socConstraints = secondordercone(Asc,b,d,gamma);
f = [-1,-2,0];
Aineq = [];
bineq = [];
Aeq = [];
beq = [];
lb = [-Inf,-Inf,0];
ub = [Inf,Inf,2];
[x,fval] = coneprog(f,socConstraints,Aineq,bineq,Aeq,beq,lb,ub,options)
Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00   0.000000e+00   5.714286e-01   1.250000e-01    0.01
   1  -7.558066e+00   0.000000e+00   7.151114e-02   1.564306e-02    0.02
   2  -7.366973e+00   0.000000e+00   1.075440e-02   2.352525e-03    0.02
   3  -8.243432e+00   0.000000e+00   5.191882e-05   1.135724e-05    0.02
   4  -8.246067e+00   0.000000e+00   2.430813e-06   5.317403e-07    0.02
   5  -8.246211e+00   0.000000e+00   6.154504e-09   1.346298e-09    0.02
Optimal solution found.
x = 3×1

    0.4851
    3.8806
    2.0000

fval = -8.2462

2 次錐計画問題の要素を作成します。時間とメモリを節約するために、最初に、最もインデックスの高い制約を作成します。

A = diag([1,2,0]);
b = zeros(3,1);
d = [0;0;1];
gamma = -1;
socConstraints(3) = secondordercone(A,b,d,gamma);
A = diag([3,0,1]);
d = [0;1;0];
socConstraints(2) = secondordercone(A,b,d,gamma);
A = diag([0;1/2;1/2]);
d = [1;0;0];
socConstraints(1) = secondordercone(A,b,d,gamma);
f = [-1;-2;-4];
options = optimoptions('coneprog','Display','iter');

problemで説明されているように、必須フィールドを含む問題構造体を作成します。

problem = struct('f',f,...
    'socConstraints',socConstraints,...
    'Aineq',[],'bineq',[],...
    'Aeq',[],'beq',[],...
    'lb',[],'ub',[],...
    'solver','coneprog',...
    'options',options);

coneprog を呼び出して問題を解きます。

[x,fval] = coneprog(problem)
Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00   0.000000e+00   5.333333e-01   5.555556e-02    0.02
   1  -9.696012e+00   5.551115e-17   7.631901e-02   7.949897e-03    0.02
   2  -1.178942e+01   0.000000e+00   1.261803e-02   1.314378e-03    0.02
   3  -1.294426e+01   9.251859e-18   1.683078e-03   1.753206e-04    0.02
   4  -1.295217e+01   9.251859e-18   8.994595e-04   9.369370e-05    0.02
   5  -1.295331e+01   0.000000e+00   4.748841e-04   4.946709e-05    0.02
   6  -1.300753e+01   0.000000e+00   2.799942e-05   2.916606e-06    0.03
   7  -1.300671e+01   9.251859e-18   2.366136e-05   2.464725e-06    0.03
   8  -1.300850e+01   9.251859e-18   8.251573e-06   8.595388e-07    0.03
   9  -1.300842e+01   4.625929e-18   7.332583e-06   7.638108e-07    0.03
  10  -1.300866e+01   9.251859e-18   2.616719e-06   2.725749e-07    0.03
  11  -1.300892e+01   1.850372e-17   2.215835e-08   2.308161e-09    0.03
Optimal solution found.
x = 3×1

    0.4238
    1.6477
    2.3225

fval = -13.0089

2 次錐計画問題を作成します。時間とメモリを節約するために、最初に、最もインデックスの高い制約を作成します。

A = diag([1,2,0]);
b = zeros(3,1);
d = [0;0;1];
gamma = -1;
socConstraints(3) = secondordercone(A,b,d,gamma);
A = diag([3,0,1]);
d = [0;1;0];
socConstraints(2) = secondordercone(A,b,d,gamma);
A = diag([0;1/2;1/2]);
d = [1;0;0];
socConstraints(1) = secondordercone(A,b,d,gamma);
f = [-1;-2;-4];
options = optimoptions('coneprog','Display','iter');
A = [];
b = [];
Aeq = [];
beq = [];
lb = [];
ub = [];

解法プロセスに関する情報を要求する問題を解きます。

[x,fval,exitflag,output] = coneprog(f,socConstraints,A,b,Aeq,beq,lb,ub,options)
Iter           Fval  Primal Infeas    Dual Infeas    Duality Gap    Time
   0   0.000000e+00   0.000000e+00   5.333333e-01   5.555556e-02    0.05
   1  -9.696012e+00   5.551115e-17   7.631901e-02   7.949897e-03    0.07
   2  -1.178942e+01   0.000000e+00   1.261803e-02   1.314378e-03    0.07
   3  -1.294426e+01   9.251859e-18   1.683078e-03   1.753206e-04    0.07
   4  -1.295217e+01   9.251859e-18   8.994595e-04   9.369370e-05    0.07
   5  -1.295331e+01   0.000000e+00   4.748841e-04   4.946709e-05    0.07
   6  -1.300753e+01   0.000000e+00   2.799942e-05   2.916606e-06    0.07
   7  -1.300671e+01   9.251859e-18   2.366136e-05   2.464725e-06    0.07
   8  -1.300850e+01   9.251859e-18   8.251573e-06   8.595388e-07    0.07
   9  -1.300842e+01   4.625929e-18   7.332583e-06   7.638108e-07    0.07
  10  -1.300866e+01   9.251859e-18   2.616719e-06   2.725749e-07    0.07
  11  -1.300892e+01   1.850372e-17   2.215835e-08   2.308161e-09    0.07
Optimal solution found.
x = 3×1

    0.4238
    1.6477
    2.3225

fval = -13.0089
exitflag = 1
output = struct with fields:
           iterations: 11
    primalfeasibility: 1.8504e-17
      dualfeasibility: 2.2158e-08
           dualitygap: 2.3082e-09
            algorithm: 'interior-point'
         linearsolver: 'augmented'
              message: 'Optimal solution found.'

  • 反復表示と出力構造体の両方に、coneprog が 12 回の反復を使用して解に到達したことが示されます。

  • 終了フラグ値 1output.message'Optimal solution found.' は、解が信頼できることを示します。

  • output 構造体には、双対性ギャップと同様に、解法プロセスを通して実行不可能性が減少する傾向にあることが示されます。

  • fval 出力は、f'*x を乗算することによって再現できます。

f'*x
ans = -13.0089

2 次錐計画問題を作成します。時間とメモリを節約するために、最初に、最もインデックスの高い制約を作成します。

A = diag([1,2,0]);
b = zeros(3,1);
d = [0;0;1];
gamma = -1;
socConstraints(3) = secondordercone(A,b,d,gamma);
A = diag([3,0,1]);
d = [0;1;0];
socConstraints(2) = secondordercone(A,b,d,gamma);
A = diag([0;1/2;1/2]);
d = [1;0;0];
socConstraints(1) = secondordercone(A,b,d,gamma);
f = [-1;-2;-4];

他のすべての coneprog 出力とともに解における双対変数を要求する問題を解きます。

[x,fval,exitflag,output,lambda] = coneprog(f,socConstraints);
Optimal solution found.

返された lambda 構造体を調査します。唯一の問題制約が錐制約のため、lambda 構造体内の soc フィールドだけを調査します。

disp(lambda.soc)
    5.9426
    4.6039
    2.4624

制約には、解において制約が有効であることを示す非ゼロの双対値が含まれています。

入力引数

すべて折りたたむ

係数ベクトル。実数ベクトルまたは実数配列として指定されます。係数ベクトルは、目的関数 f'*x を表します。表記では、f が列ベクトルになっていますが、行ベクトルや配列も使用できます。coneprog は配列 f を列ベクトル f(:) に内部的に変換します。

例: f = [1,3,5,-6]

データ型: double

SecondOrderConeConstraint オブジェクトのベクトルとして指定された、2 次錐制約secondordercone 関数を使用してこれらのオブジェクトを作成します。

socConstraints が次の制約をエンコードします。

Asc(i)xbsc(i)dscT(i)xγ(i)

ここで、配列と方程式の対応付けを以下に示します。

  • Asc(i) = socConstraints.A(i)

  • bsc(i) = socConstraints.b(i)

  • dsc(i) = socConstraints.d(i)

  • γ(i) = socConstraints.gamma(i)

例: Asc = diag([1 1/2 0]); bsc = zeros(3,1); dsc = [0;0;1]; gamma = -1; socConstraints = secondordercone(Asc,bsc,dsc,gamma);

データ型: struct

実数行列として指定される線形不等式制約です。AMN 列の行列で、M は不等式の数、N は変数の数 (f の長さ) です。大規模な問題の場合は、A をスパース行列として渡します。

AM 個の線形不等式を符号化します。

A*x <= b,

ここで、xN 個の変数 x(:) の列ベクトル、bM 個の要素をもつ列ベクトルです。

たとえば、次の不等式を考えてみましょう。

x1 + 2x2 ≤ 10
3x1 + 4x2 ≤ 20
5x1 + 6x2 ≤ 30.

次の制約を入力することによって、不等式を指定します。

A = [1,2;3,4;5,6];
b = [10;20;30];

例: x 成分の和が 1 以下になるように指定するために、A = ones(1,N)b = 1 をとります。

データ型: double

実数ベクトルで指定される線形不等式制約です。b は、行列 A に関連する M 要素ベクトルです。b を行ベクトルとして渡す場合、ソルバーは b を列ベクトル b(:) に内部的に変換します。大規模な問題の場合は、b をスパース ベクトルとして渡します。

bM 個の線形不等式を符号化します。

A*x <= b,

ここで、xN 個の変数 x(:) の列ベクトル、AMN 列の行列です。

たとえば、次の不等式を考えてみましょう。

x1 + 2x2 ≤ 10
3x1 + 4x2 ≤ 20
5x1 + 6x2 ≤ 30.

次の制約を入力することによって、不等式を指定します。

A = [1,2;3,4;5,6];
b = [10;20;30];

例: x の成分の和が 1 以下であることを指定するには、A = ones(1,N)b = 1 を使用します。

データ型: double

実数行列として指定される線形等式制約です。AeqMeN 列の行列で、Me は等式の数、N は変数の数 (f の長さ) です。大規模な問題の場合は、Aeq をスパース行列として渡します。

AeqMe 個の線形等式を符号化します。

Aeq*x = beq,

ここで、xN 個の変数 x(:) の列ベクトル、beqMe 個の要素をもつ列ベクトルです。

たとえば、次の等式を考えてみましょう。

x1 + 2x2 + 3x3 = 10
2x1 + 4x2 + x3 = 20.

次の制約を入力することによって、等式を指定します。

Aeq = [1,2,3;2,4,1];
beq = [10;20];

例: x 成分の和が 1 になるように指定するために、Aeq = ones(1,N)beq = 1 をとります。

データ型: double

実数ベクトルで指定される線形等式制約です。beq は、行列 Aeq に関連する Me 要素ベクトルです。beq を行ベクトルとして渡す場合、ソルバーは beq を列ベクトル beq(:) に内部的に変換します。大規模な問題の場合は、beq をスパース ベクトルとして渡します。

beqMe 個の線形等式を符号化します。

Aeq*x = beq,

ここで、xN 個の変数 x(:) の列ベクトル、AeqMeN 列の行列です。

たとえば、次の等式を考えてみましょう。

x1 + 2x2 + 3x3 = 10
2x1 + 4x2 + x3 = 20.

次の制約を入力することによって、等式を指定します。

Aeq = [1,2,3;2,4,1];
beq = [10;20];

例: x の成分の和が 1 であることを指定するには、Aeq = ones(1,N)beq = 1 を使用します。

データ型: double

下限。実数ベクトルまたは実数配列として指定されます。f の長さが lb の長さと等しい場合、lb は次を指定します。

x(i) >= lb(i) (すべての i について)

numel(lb) < numel(f) の場合、lb は次を指定します。

x(i) >= lb(i) (1 <= i <= numel(lb))

この場合、ソルバーによって警告が発行されます。

例: すべての x 成分が正になるように指定するには、lb = zeros(size(f)) を使用します。

データ型: double

実数ベクトルまたは実数配列として指定される上限です。f の長さが ub の長さと等しい場合、ub は次を指定します。

x(i) <= ub(i) (すべての i について)

numel(ub) < numel(f) の場合、ub は次を指定します。

x(i) <= ub(i) (1 <= i <= numel(ub))

この場合、ソルバーによって警告が発行されます。

例: すべての x 成分が 1 未満になるように指定するには、ub = ones(size(f)) を使用します。

データ型: double

最適化オプション。optimoptions の出力として指定されます。

オプション説明
ConstraintTolerance

制約の実行可能性の許容誤差で 0 から 1 のスカラー。ConstraintTolerance は、主実行可能性の許容誤差を測定します。既定値は 1e-6 です。

Display

表示レベル (反復表示を参照):

  • 'final' (既定) — 最終出力のみを表示する。

  • 'iter' — 各反復の出力を表示する。

  • 'off' または 'none' — 出力を表示しない。

LinearSolver

反復における 1 ステップを解くためのアルゴリズム:

  • 'auto' (既定) — coneprog がステップ ソルバーを選択する。

    • スパースな問題の場合は 'prodchol' をステップ ソルバーとする。

    • それ以外の場合は 'augmented' をステップ ソルバーとする。

  • 'augmented' — 拡張形式のステップ ソルバー。詳細については、[1]を参照してください。

  • 'normal' — 標準形式のステップ ソルバー。詳細については、[1]を参照してください。

  • 'normal-dense' — 密な線形代数を使用した標準形式のステップ ソルバー。

  • 'prodchol' — 積形式のコレスキー ステップ ソルバー。詳細については、[4][5]を参照してください。

  • 'schur' — シューア補行列法ステップ ソルバー。詳細については、[2]を参照してください。

'auto' がうまく動作しなかった場合は、LinearSolver で次の方法を試してみてください。

  • スパースな問題の場合は、'normal' を試してみます。

  • スパースな問題で、なおかつ密な列がいくつかある場合や大きな円錐がある場合は、'prodchol' または 'schur' を試してみます。

  • 密な問題の場合は、'normal-dense' または 'augmented' を使用します。

スパースの例については、coneprog アルゴリズムの速度の比較を参照してください。

MaxIterations

反復の最大許容回数 (非負の整数)。既定値は 200 です。

詳細については、許容誤差と停止条件反復と関数カウントを参照してください。

MaxTime

アルゴリズムの秒単位の最大実行時間 (非負の数値または Inf)。既定値は Inf で、この停止条件を無効にします。

OptimalityTolerance

双対実行可能性に関する終了許容誤差 (非負のスカラー)。既定値は 1e-6 です。

例: optimoptions('coneprog','Display','iter','MaxIterations',100)

次のフィールドをもつ構造体として指定される問題構造体です。

フィールド名エントリ

f

線形目的関数ベクトル f

socConstraints

2 次錐制約の構造体配列

Aineq

線形不等式制約の行列

bineq

線形不等式制約のベクトル

Aeq

線形等式制約の行列

beq

線形等式制約のベクトル
lb下限のベクトル
ub上限のベクトル

solver

'coneprog'

options

optimoptions で作成されたオプション

データ型: struct

出力引数

すべて折りたたむ

実数ベクトルまたは実数配列として返される解です。x のサイズは、f のサイズと同じです。x の出力は、exitflag の値が -2、-3、または -10 の場合に空になります。

解での目的関数値。実数として返されます。一般的に、fval = f'*x になります。fval の出力は、exitflag の値が -2、-3、または -10 の場合に空になります。

coneprog の停止理由。整数として返されます。

説明

1

関数が解 x に収束したことを示します。

0

反復回数が options.MaxIterations を超過したか、秒単位の求解時間が options.MaxTime を超過しました。

-2

実行可能な点が見つかりません。

-3

問題が非有界です。

-7

探索方向が小さくなりすぎ、これ以上進むことができないことを示します。

-10

問題が数値的に不安定です。

ヒント

終了フラグ値として 0-7、または -10 を取得した場合は、LinearSolver オプションの値を変えてみてください。

最適化プロセスに関する情報。次のフィールドをもつ構造体として返されます。

フィールド説明
algorithm

使用される最適化アルゴリズム

dualfeasibility

双対制約違反の最大値

dualitygap

双対性ギャップ

iterations

反復回数

message

終了メッセージ

primalfeasibility

制約違反の最大値

linearsolver使用される内部ステップ ソルバー アルゴリズム

output のフィールドの dualfeasibilitydualitygap、および primalfeasibility は、exitflag の値が –2、–3、または –10 の場合に空になります。

解における双対変数。次のフィールドをもつ構造体として返されます。

フィールド説明
lower

lb に対応する下限

upper

ub に対応する上限

ineqlin

A および b に対応する線形不等式

eqlin

Aeq および beq に対応する線形等式

socsocConstraints に対応する 2 次錐制約

lambda は、exitflag の値が -2、-3、または -10 の場合に空 ([]) になります。

ラグランジュ乗数 (双対変数) は、解において一定 (ゼロ勾配) の次のラグランジュ関数の一部です。

fTx+iλsoc(i)(Asoc(i)xbsoc(i)dsocT(i)x+gamma(i))      +λineqlinT(Axb)+λeqlinT(Aeqxbeq)+λupperT(xub)+λlowerT(lbx).

詳細

すべて折りたたむ

2 次錐制約

制約の

AxbdTxγ

が 2 次錐制約と呼ばれる理由3 次元空間で x-y 平面が楕円形の断面で、直径が z 座標に直交している円錐を考えます。y 座標のスケールは ½ で、x座標のスケールは 1 です。[0,0,0] にある点でこの円錐の内部を定義する不等式は次のとおりです。

x2+y24z.

coneprog 構文では、この円錐は次の引数をとります。

A = diag([1 1/2 0]);
b = [0;0;0];
d = [0;0;1];
gamma = 0;

円錐の境界をプロットします。

[X,Y] = meshgrid(-2:0.1:2);
Z = sqrt(X.^2 + Y.^2/4);
surf(X,Y,Z)
view(8,2)
xlabel 'x'
ylabel 'y'
zlabel 'z'

Cone with point at zero, widening in the vertical direction.

b 引数と gamma引数によって、円錐が移動します。A 引数と d引数によって、円錐が回転し、その形状が変化します。

アルゴリズム

アルゴリズムは内点法を使用します。詳細については、2 次錐計画法アルゴリズムを参照してください。

代替機能

アプリ

[最適化] ライブ エディター タスクが coneprog にビジュアル インターフェイスを提供します。

バージョン履歴

R2020b で導入

すべて展開する