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

二次制約付きの線形または二次目的関数

この例では、線形または二次の目的関数および二次不等式制約をもつ最適化問題を解く方法を説明します。ここでは、目的関数と制約関数の勾配およびヘッシアンを生成して使用する方法を説明します。

二次制約付きの問題

問題を次のように表現できると仮定します。

minx(12xTQx+fTx+c)

制約条件

12xTHix+kiTx+di0,

ここで、1 ≤ i ≤ m です。少なくとも 1 つの Hi が非ゼロであると仮定します。そうでない場合は quadprog または linprog を使用してこの問題を解くことができます。Hi が非ゼロの場合は制約が非線形になり、最適化の意思決定表に示すように、適切なソルバーは fmincon です。

この例では二次行列を対称行列と仮定しています。これによる一般性損失はありません。非対称行列 H (または Q) は、同等な対称化された行列 (H + HT)/2 による置き換えが可能です。

x に N 個の成分がある場合、Q と Hi は N 行 N 列の行列、f と ki は N 行 1 列のベクトル、c と di はスカラーになります。

目的関数

fmincon 構文を使用して問題を定式化します。xf は列ベクトルであると仮定します (初期ベクトル x0 が列ベクトルの場合、x は列ベクトルになります)。

function [y,grady] = quadobj(x,Q,f,c)
y = 1/2*x'*Q*x + f'*x + c;
if nargout > 1
    grady = Q*x + f;
end

制約関数

一貫性を保ち、インデックス処理を容易にするため、各二次制約を 1 つの cell 配列に配置します。同様に、線形項および定数項をそれぞれ cell 配列に配置します。

function [y,yeq,grady,gradyeq] = quadconstr(x,H,k,d)
jj = length(H); % jj is the number of inequality constraints
y = zeros(1,jj);
for i = 1:jj
    y(i) = 1/2*x'*H{i}*x + k{i}'*x + d{i};
end
yeq = [];
    
if nargout > 2
    grady = zeros(length(x),jj);
    for i = 1:jj
        grady(:,i) = H{i}*x + k{i};
    end
end
gradyeq = [];

数値の例

たとえば、以下の問題について考えてみましょう。

Q = [3,2,1;
     2,4,0;
     1,0,5];
f = [-24;-48;-130];
c = -2;

rng default % for reproducibility
% Two sets of random quadratic constraints:
H{1} = gallery('randcorr',3); % random positive definite matrix
H{2} = gallery('randcorr',3);
k{1} = randn(3,1);
k{2} = randn(3,1);
d{1} = randn;
d{2} = randn;

ヘッシアン

ヘッセ関数を作成します。ラグランジュ関数のヘッシアンは次の式で与えられます。

xx2L(x,λ)=2f(x)+λi2ci(x)+λi2ceqi(x).

fmincon は、ラグランジュ乗数 λi の近似セットを計算し、これらを構造体にパッケージ化します。ヘッシアンを含めるには、次の関数を使用します。

function hess = quadhess(x,lambda,Q,H)
hess = Q;
jj = length(H); % jj is the number of inequality constraints
for i = 1:jj
    hess = hess + lambda.ineqnonlin(i)*H{i};
end

問題を最も効率的に解くには、fmincon interior-point アルゴリズムを使用します。このアルゴリズムは、ユーザーが指定したヘッセ関数を受け入れます。以下のオプションを設定します。

options = optimoptions(@fmincon,'Algorithm','interior-point',...
    'SpecifyObjectiveGradient',true,'SpecifyConstraintGradient',true,...
    'HessianFcn',@(x,lambda)quadhess(x,lambda,Q,H));

fmincon を呼び出してこの問題を解きます。

fun = @(x)quadobj(x,Q,f,c);
nonlconstr = @(x)quadconstr(x,H,k,d);
x0 = [0;0;0]; % column vector
[x,fval,eflag,output,lambda] = fmincon(fun,x0,...
    [],[],[],[],[],[],nonlconstr,options);

ラグランジュ乗数を確認します。

lambda.ineqnonlin
ans =

   12.8412
   39.2337

非線形不等式の乗数がどちらも非ゼロなので、この解では両方の二次制約がアクティブです。

ヘッシアンを使用する際の効率性

勾配とヘッシアンを使った内点法アルゴリズムは効率的です。関数評価の回数を調べてみましょう。

output
output = 

         iterations: 9
          funcCount: 10
    constrviolation: 0
           stepsize: 5.3547e-04
          algorithm: 'interior-point'
      firstorderopt: 1.5851e-05
       cgiterations: 0
            message: 'Local minimum found that satisfies the constraints.

Optimization compl...'

fmincon は関数評価を 10 回実行しただけで問題を解いています。

これを、ヘッシアンを使用せずに問題を解く場合と比較してみましょう。

options.HessianFcn = [];
[x2,fval2,eflag2,output2,lambda2] = fmincon(fun,[0;0;0],...
    [],[],[],[],[],[],nonlconstr,options);
output2
output2 = 

         iterations: 17
          funcCount: 22
    constrviolation: 0
           stepsize: 2.8475e-04
          algorithm: 'interior-point'
      firstorderopt: 1.7680e-05
       cgiterations: 0
            message: 'Local minimum found that satisfies the constraints.

Optimization compl...'

この場合、fmincon が約 2 倍の回数の反復処理と関数評価を実行しています。解は許容誤差範囲内において同一です。

二次等式制約の拡張

二次等式制約もある場合、基本的に上記と同様の手法を使用できます。問題は同じですが、制約が追加されます。

12xTJix+piTx+qi=0.

Ji、pi、qi の各変数を使用するように、制約を再定式化します。lambda.eqnonlin(i) 構造体に等式制約のためのラグランジュ乗数が含まれています。

関連する例

詳細