Main Content

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

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

二次制約付きの問題

問題は次の形式をとるものとします。

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) 構造体に等式制約のためのラグランジュ乗数が含まれています。

関連するトピック