最適化パッケーコマンドにおける行列の入力
|
行列の形の使用
|
|
•
|
およびが行列およびベクトルの形で与えられる問題、もしくはベクトルおよび行列の手続きで与えられる問題のとき、ソルバは最も効率よく働きます。行列の形は前処理が最も少ないことが必要とされます。なぜなら、内部のソルバによって使われる入力がほとんど似ているためです。前処理を減らし早く計算することができます。
|
|
最小化 f(x)
制約条件
v(x) <= 0 (非線形不等式制約条件)
w(x) = 0 (非線形等式制約条件)
A x <= b (線形不等式制約条件)
Aeq x = beq (線形等式制約条件)
bl <= x <= bu (境界条件)
|
|
ただし、x は変数ベクトル;f(x) は実目的関数;v(x) と w(x) は x の関数ベクトル;b, beq, bl および bu はベクトル;A および Aeq は行列です。x の次数は n です。行列とベクトルは要素ごとに対応しています。各要素の詳細は次の章で述べます。
|
|
|
目的関数
|
|
|
線形目的関数
|
|
•
|
LPSolve コマンドは線形目的関数 f(x) を係数を要素とする n 次元ベクトル c で指定することができます。
|
|
|
2 次目的関数
|
|
•
|
QPSolve コマンドは 2 次目的関数を扱うことができます。目的関数は c をベクトル、H を対称行列とするとき、c' x + (1/2) x' H x の形であることを仮定します。線形関数が含まれている場合は [c, H] と、そうでない場合は、H と指定します。H は n 次正方行列、c は n 次元ベクトルである必要があります。 c' および x' はそれぞれ c と x の転置ベクトルをあらわします。
|
|
|
非線形目的関数
|
|
•
|
NLPSolve コマンドは procedure[float](Vector) タイプの手続き p で作られた非線形目的関数 f(x) を使うことができます。この手続きは x であらわされる n 次のパラメータベクトルを入力し、f(x) の値を返します。n の値は手続き p から決めることができないため NLPSolve コマンドに入力します。
|
•
|
計算を最適にするには、objectivegradient=objgrd を使って目的関数のを入力します。ただし、objgrd は型 procedure(Vector, Vector) の手続きで作成します。この手続きは x であらわされるパラメータベクトルを入力し、x による f の勾配を計算し第 2 引数にベクトルパラメータを返します。各ベクトルの次元は n です。
|
|
|
線形最小二乗目的関数
|
|
•
|
LSSolve コマンドは線形最小二乗目的関数を使うことができます。f(x) を (1/2) ||c - G x||^2 で与えられる形最小二乗目的関数とします。このとき f(x) は c を q 次元ベクトル、G を q × n 行列とするとき [c, G] のリストで指定します。
|
|
|
非線形最小二乗目的関数
|
|
•
|
LSSolve コマンドは非線形最小二乗目的関数を使うことができます。f(x) は (1/2) (f1(x)^2 + f2(x)^2 + ... + fq(x)^2) で与えられたとします。f(x) は型 procedure(Vector, Vector) を持つ手続き p で与えます。この手続きは x であわらされる n 次元パラメータベクトルを入力し、x による残差 f1, f2, ..., fq を計算し、第 2 引数のベクトルパラメータにこれらの残差を返します。n と q の値は、手続き p から決定できないので LSSolve コマンドに個々に入力します。
|
•
|
手続き p は第 3 引数にタイプ float を持つ needfi を指定することができます。この場合、p は残差の計算が終わるまで needfi < 0.0 という制約を受けます、もしくは i が trunc(needfi) となる i 番目の残差のみが計算されます 。これにより必要な残差のみ計算が行われ効率よく計算できます。このように手続きは次のような形です:
|
•
|
p := proc(x, y, needfi)
if needfi < 0.0 or trunc(needfi) = 1 then
y[1] := ...;
end if;
if needfi < 0.0 or trunc(needfi) = 2 then
y[2] := ...;
end if;
...
end proc;
|
•
|
計算を最適に行うには、 objectivejacobian=objjac オプションを使って目的関数の を入力します。ただし、objjac は procedure(Vector, Matrix) のタイプを持ちます。この手続きは x に関するパラメータベクトルを入力し x に関する f のヤコビ行列計算し、 q × n パラメータ行列を返します。
|
|
|
|
制約条件
|
|
|
線形制約条件
|
|
•
|
線形制約条件は行列の入力を扱うすべての Optimization ソルバで使うことができます。次の形の中から使うものを選びます。
|
|
[A, b, Aeq, beq] -- 不等式および等式制約条件
|
|
[NoUserValue, NoUserValue, Aeq, beq] -- 等式制約条件のみ
|
•
|
A および b は線形不等式制約条件A x <= b の要素をあらわし、Aeq および beq は線形等式制約条件 Aeq x = b の要素をあらわします。
|
•
|
b および beq はベクトルまたは数値のタイプを持つことができます。制約条件が 1 つの数値ならば、要素すべてがその数値であるベクトルを表現します。
|
•
|
A および Aeq は列が n の行列でなければなりません。A の行の長さはベクトル b と同じである必要があり、Aeq の行の長さはベクトル beq と同じである必要があります。
|
|
|
非線形制約条件
|
|
•
|
非線形制約条件は NLPSolve および LSSolve コマンドで使うことができます。非線形不等式制約条件 v(x) <= 0 および非線形等式制約条件 w(x) = 0 を考えます。y(x) を w(x) に従う v(x) からなるベクトルとする。y(x) は型 procedure(Vector, Vector) の手続き pcons にあらわされます。ただし、第 1 引数は現在の点 x 、第 2 引数は点 x での y の値をあらわすパラメータです。Optimization に含まれるコマンドは制約式の数を数えることはできないので、数は使い方のように指定しなければいけません。
|
•
|
手続き pcons は制約条件の数を次元とする数値ベクトルを第 3 引数 needc を持つことができます。この場合、p は needc[i] の値でチェックを行い、 0.0 より大きければ i 番目の制約条件の左辺を計算します。これにより、必要なときだけ制約条件を使うので計算は効率的になります。
|
•
|
pcons := proc(x, y, needc)
if needc[1] > 0.0 then
y[1] := ...;
end if;
if needc[2] > 0.0 then
y[2] := ...;
end if;
...
end proc;
|
•
|
計算を最適にするためには、constraintjacobian=conjac オプションを使い y(x) のを指定します。ただし、conjac は型 procedure(Vector, Matrix) を持ちます。この手続きは x によるベクトルによって x による y のヤコビ行列を計算し、行列を値として返します。上で説明した手続き pcons のように手続き conjac は第 3 引数として needc を使うことができます。needc は needc[i]>0 のとき i 番目のヤコビ行列を評価するような数値ベクトルです。
|
|
|
|
境界条件
|
|
•
|
変数の境界はすべての Optimization ソルバで行列の形の入力をすることができます。境界とは上界および下界からなるリスト [bl, bu] です。もしリストが空ならば、境界条件は実数空間となります。
|
•
|
境界 bl および bu は n 次ベクトルまたは数値で指定することができます。数値が指定されたときは、要素すべてがその数値である n 次ベクトルを表現します。
|
•
|
上界はあるが下界がない場合、bl は -Float(infinity) になります。下界はあるが上界がない場合、bu は Float(infinity) になります。
|
|
|
初期値
|
|
•
|
初期値は行列を入力として許すすべての Optimization ソルバで initialpoint=p オプションによって使うことができます。ただし、p は n 次ベクトルです。initialpoint オプションについて詳しくは Optimization/Options を参照ください。
|
|
|
解
|
|
•
|
objval を最小 (最大) 値とする浮動小数、 solpt をその点 () に対応するベクトルとするとき、解はリスト [objval, solpt] で与えられます。
|
|
|
Additional Requirements追加条件
|
|
•
|
最適な計算をするためには、要素が次の追加条件を満たすことが必要です。
行列およびベクトルは datatype=float を指定して構成します。そうでないとき、ソルバは変換を行わなければいけません。
QPSolve コマンドの行列のような対称行列は、 shape=symmetric および storage=rectangular を使い構成すべきです。そうでないとき、ソルバは変換を行わなければいけません。
入力手続きは evalhf コマンドによって評価したもので構成するべきです。そうでないとき、ソルバはより効率的な evalhf での計算を利用することはできません。
スカラーを返す手続きはすべての場合において表示は簡単になりませんが、浮動小数で結果を返すべきです。
|
•
|
次のは必ずしも必要ではありまえん。出力パラメータを使いベクトルや行列を返す手続きであるとき、パラメータに対して浮動小数の値を割り当てることができます。そうでない場合はエラーを起こします。
|
|
|
例
|
|
行列の形で線形計画問題を記述し、Optimization[LPSolve] を使って解きます。
>
|
c := Vector([-1,-1], datatype=float):
A := Matrix([[-3,1],[5,1]], datatype=float):
b := Vector([0.5,2], datatype=float):
LPSolve(c, [A, b], [0.0, infinity]);
|
| (8.1) |
行列の形で 2 次計画問題を記述し、Optimization[QPSolve] を使って解きます。
>
|
c := Vector([2,5], datatype=float):
H := Matrix([[6,3],[3,4]], datatype=float):
A := Matrix([[-1,1]], datatype=float):
b := Vector([-2], datatype=float):
QPSolve([c, H], [A, b]);
|
| (8.2) |
非線形制約条件 w^2+x^2+y+z<=5 および線形制約条件 z=v のもと非線形目的関数 w^3*(v-w)^2+(w-x-1)^2+(x-y-2)^2+(y-z-3)^2 について考えます。 V は変数ベクトル <v,x,w,y,z> とするとき、 V を用いて目的関数を浮動小数の値をとる手続きとして書きます。
>
|
p := proc (V)
V[3]^3*(V[1]-V[3])^2+(V[3]-V[2]-1)^2+
(V[2]-V[4]-2)^2+(V[4]-V[5]-3)^2
end proc:
|
非線形制約条件 w^2+x^2+y+z<=5 を入力パラメータ V と出力パラメータとしてベクトル W を用いて手続きを書きます。needc をパラメータとして入力します。
>
|
nlc := proc(V, W, needc)
if needc[1] > 0.0 then
W[1] := V[3]^2+V[2]^2+V[4]+V[5]-5
end if:
end proc:
|
線形制約条件をあらわすリストを構成します。この場合、制約条件は z-u=0 です。
>
|
Aeq := Matrix([[-1,0,0,0,1]], datatype=float):
beq := Vector([0], datatype=float):
lc := [NoUserValue, NoUserValue, Aeq, beq]:
|
よりよい計算のためには、目的関数の勾配および制約条件のヤコビ行列をあらわす手続きを使います。手続き objgra は V を入力パラメータとしてベクトルパラメータ W にV の勾配を返します。
>
|
objgra := proc (V, W)
W[1] := 2*V[3]^3*(V[1]-V[3]);
W[2] :=-2*V[3]+4*V[2]-2-2*V[4];
W[3] := 3*V[3]^2*(V[1]-V[3])^2-2*V[3]^3*(V[1]-V[3])+
2*V[3]-2*V[2]-2;
W[4] :=-2*V[2]+4*V[4]-2-2*V[5];
W[5] :=-2*V[4]+2*V[5]+6
end proc:
|
手続き conjac は V を入力パラメータとして出力行列パラメータ M にヤコビ行列を返します。
>
|
conjac := proc(V, M, needc)
if needc[1] > 0.0 then
M[1,1] := 0;
M[1,2] := 2*V[2];
M[1,3] := 2*V[3];
M[1,4] := 1;
M[1,5] := 1;
end if
end proc:
|
Optimization[NLPSolve] コマンドの第 1 引数に変数の数 5 を入力し非線形計画問題を解きます。変数に非負条件を指定するため assume=nonnegative オプションを使います。
>
|
NLPSolve(5, p, 1, nlc, lc, objectivegradient=objgra, constraintjacobian=conjac, assume=nonnegative);
|
| (8.3) |
|
|