制約をもつ非線形システム

不等式制約をもつ方程式の解法

fsolve は非線形方程式系を解きます。しかし、範囲も含めて、制約を一切使用できません。では、制約が含まれる非線形方程式系を解くにはどうすればよいでしょうか。

制約を満たす解が存在するという保証はありません。実際、制約を満たさない解も含め、問題に解が存在しない可能性もあります。しかしながら、制約を満たす解を求めるために役立つ手法があります。

この手法を説明するために、次の方程式の解き方を考えてみましょう。

$\begin{array}{c} F_{1} (x) = (x_{1} + 1) (10 - x_{1}) \frac{1 + x_{2}^{2}}{1 + x_{2}^{2} + x_{2}} \\ F_{2} (x) = (x_{2} + 2) (20 - x_{2}) \frac{1 + x_{1}^{2}}{1 + x_{1}^{2} + x_{1}}, \end{array}$

ここで、 $x$ の成分は非負でなければなりません。これらの方程式には次の 4 つの解があります。

$\begin{array}{l} x = (- 1, - 2) \\ x = (10, - 2) \\ x = (- 1, 20) \\ x = (10, 20) . \end{array}$

制約を満たす解は $x = (10, 20)$ の 1 つだけです。

この例の最後に記載されている fbnd 補助関数が、 $F (x)$ を数値的に計算します。

異なる開始点の使用

一般に、 $N$ 変数の $N$ 次方程式系では解が孤立しています。つまり、各解の近傍に他の解はありません。このため、一部の制約を満たす解の探究方法の 1 つは、いくつかの初期点 x0 を生成し、それぞれの x0 から fsolve の実行を開始することです。

この例では、方程式系 $F (x) = 0$ の解を求めるために、平均 0、標準偏差 100 で正規分布する 10 個の無作為な点を使用します。

rng default % For reproducibility
N = 10; % Try 10 random start points
pts = 100*randn(N,2); % Initial points are rows in pts
soln = zeros(N,2); % Allocate solution
opts = optimoptions('fsolve','Display','off');
for k = 1:N
    soln(k,:) = fsolve(@fbnd,pts(k,:),opts); % Find solutions
end

制約を満たす解を列挙します。

idx = soln(:,1) >= 0 & soln(:,2) >= 0;
disp(soln(idx,:))

   10.0000   20.0000
   10.0000   20.0000
   10.0000   20.0000
   10.0000   20.0000
   10.0000   20.0000

異なるアルゴリズムの使用

fsolve は次の 3 つのアルゴリズムをもっています。その各々で異なる解が見つかります。

この例では、x0 = [1,9] を使用して各アルゴリズムが返す解を調べます。

x0 = [1,9];
opts = optimoptions(@fsolve,'Display','off',...
    'Algorithm','trust-region-dogleg');
x1 = fsolve(@fbnd,x0,opts)

x1 = 1×2

   -1.0000   -2.0000

opts.Algorithm = 'trust-region';
x2 = fsolve(@fbnd,x0,opts)

x2 = 1×2

   -1.0000   20.0000

opts.Algorithm = 'levenberg-marquardt';
x3 = fsolve(@fbnd,x0,opts)

x3 = 1×2

    0.9523    8.9941

ここでは、同じ初期点に対し 3 つのアルゴリズムすべてで異なる解が見つかりました。いずれの解も制約を満たしていません。報告された "解" x3 は、解でさえなく、局所的な停留点にすぎません。

範囲を指定した `lsqnonlin` の使用

lsqnonlin は、ベクトル関数 $F (x)$ の成分の二乗和を最小化しようとします。このため、方程式 $F (x) = 0$ の解を求めようと試みます。また、lsqnonlin は範囲の制約を受け入れます。

lsqnonlin 用に例の問題を定式化して解を求めます。

lb = [0,0];
rng default
x0 = 100*randn(2,1);
[x,res] = lsqnonlin(@fbnd,x0,lb)

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

res = 2.4783e-25

この場合は、lsqnonlin が制約を満たす解に収束します。lsqnonlin を Global Optimization Toolbox の MultiStart ソルバーと併用して、多くの初期点を自動的に調べることができます。詳細は、MultiStart Using lsqcurvefit or lsqnonlin (Global Optimization Toolbox)を参照してください。

`fmincon` の制約としての方程式と不等式の設定

次のように、問題を再定式化して fmincon を使用できます。

すべての x に対して 0 に評価される @(x)0 などの定数目的関数を指定します。
fmincon で fsolve 目的関数を非線形等式制約として設定します。
通常の fmincon 構文で他の制約を指定します。

この例の最後にある fminconstr 補助関数が非線形制約を実装します。制約付きの問題を解きます。

lb = [0,0]; % Lower bound constraint
rng default % Reproducible initial point
x0 = 100*randn(2,1);
opts = optimoptions(@fmincon,'Algorithm','interior-point','Display','off');
x = fmincon(@(x)0,x0,[],[],[],[],lb,[],@fminconstr,opts)

この場合は、fmincon が開始点から問題を解きます。

補助関数

次のコードは、補助関数 fbnd を作成します。

function F = fbnd(x)

F(1) = (x(1)+1)*(10-x(1))*(1+x(2)^2)/(1+x(2)^2+x(2));
F(2) = (x(2)+2)*(20-x(2))*(1+x(1)^2)/(1+x(1)^2+x(1));
end

次のコードは、補助関数 fminconstr を作成します。

function [c,ceq] = fminconstr(x)

c = []; % No nonlinear inequality
ceq = fbnd(x); % fsolve objective is fmincon nonlinear equality constraints
end

参考

fsolve | lsqnonlin | fmincon