最適化計算のためのWolfram言語12の新機能

丸山耕司
Wolfram Research Inc.

v12での新機能

凸最適化の最新ソルバー

LinearOptimization — minimize {c.x|a.x+b \[VectorGreaterEqual]0}

LinearFractionalOptimization — minimize {(α.x+β)/(γ.x+δ|a.x+b\[VectorGreaterEqual]0}

QuadraticOptimization — minimize {(1/2)x.q.x+c.x|a.x+b\[VectorGreaterEqual]0}

SecondOrderConeOptimization — minimize {c.x|‖(Subscript[a, j]).x+(Subscript[b, j])‖≤(Subscript[α, j])+(Subscript[β, j]),j=1… k}

SemidefiniteOptimization — minimize {c.x|(Subscript[a, 0])+(Subscript[x, 1])(Subscript[a, 1])+…+(Subscript[x, k])(Subscript[a, k])(Subscript[\[VectorGreaterEqual], (Subscript[𝒮, +]^n)])0}

ConicOptimization — minimize {c.x|(Subscript[a, j]).x+(Subscript[b, j])(Subscript[\[VectorGreaterEqual], (Subscript[κ, j])])0, j=1… k}

... これまでもあった?

汎用最適化ソルバー

FindMinimum —制約条件つき局所的数値最適化

NMinimize — 制約条件つき大域的数値最適化

これらは新しくない．
Wolfram言語12では，問題が凸であると(シンボリックに)認識され，ソルバーが理解可能な形式に変換されるとき，これまでより効率的に(わかりやすい形で)求解される．

フィッティング

FitRegularization — Fit and FindFitのための「正規化」, e.g., LASSO, Variation, Curvature, etc.

v12での新機能

モデリングツール

ベクトル不等式

次元の推定を伴うベクトル変数の利用

Parameter “constraint” equations

効率的な標準形への変換

ベクトル不等式

ベクトルの成分ごとの比較がすべて「真」であれば，真(True)を返す．そうでなければ，すべて偽(False)．
(よって，ベクトル空間のなかに半順序を定める)

VectorLess  {(Subscript[x, 1]),…,(Subscript[x, n])}\[VectorLess]{(Subscript[y, 1]),…,(Subscript[y, n])} if  (Subscript[x, i])<(Subscript[y, i]) for  i=1 … n

VectorLessEqual  {(Subscript[x, 1]),…,(Subscript[x, n])}\[VectorLessEqual]{(Subscript[y, 1]),…,(Subscript[y, n])} if  (Subscript[x, i])≤(Subscript[y, i]) for  i=1 … n

VectorGreaterEqual  {(Subscript[x, 1]),…,(Subscript[x, n])}\[VectorGreaterEqual]{(Subscript[y, 1]),…,(Subscript[y, n])} if  (Subscript[x, i])≥(Subscript[y, i]) for  i=1 … n

VectorGreater  {(Subscript[x, 1]),…,(Subscript[x, n])}\[VectorGreater]{(Subscript[y, 1]),…,(Subscript[y, n])} if  (Subscript[x, i])>(Subscript[y, i]) for  i=1 … n

Manipulate[Quiet@RegionPlot[{x,y}\[VectorGreaterEqual]p,{x,-2,2},{y,-2,2}],{{p,{0,0}},Locator}]

{1,2}\[VectorLess]{3,4}

どちらかがスカラーでもかまわない． (矛盾なく比較できるならどんな次元でもOK):

{1,2}\[VectorLess]a

{1,2}\[VectorLess]a /. a->3

{1,2}\[VectorLess]a /. a->{{1,2,3},{4,5}}

一般化されたベクトル不等式

ふたつのベクトルx, yについて，x-y∈κ が成り立つなら不等式 x(Underscript[\[VectorGreaterEqual], κ])y が真．κ は凸錐(convex cone)．

“NonNegativeCone” ((NonNegativeConeList[n])) x∈(ℝ^n) such that (Subscript[x, i])≥0 (“第１象限”）

“NormCone” ((NormConeList[n])) x∈(ℝ^n) such that Norm[{(Subscript[x, 1]),…,(Subscript[x, n-1])}]≤(Subscript[x, n]) (2次錐)

“SemidefiniteCone” ((SemidefiniteConeList[n]))  symmetric positive semidefinite matrices x∈(ℝ^(n×n)) (対称半正定値行列の錐)

錐 (cone)

ベクトル空間Vの部分空間Cをとり, Cの要素xをα倍 (α > 0) したαxもCの要素であるようなとき，Cは錐であるという．

2次錐不等式

“NormCone” ((NormConeList[n])) x∈(ℝ^n) such that Norm[{(Subscript[x, 1]),…,(Subscript[x, n-1])}]≤(Subscript[x, n])

Quiet[RegionPlot[{x,y,1}(Subscript[\[VectorGreaterEqual], {"NormCone",3}])0,{x,-1,1},{y,-1,1}]]

RegionPlot3D[{x,y,z}(Underscript[\[VectorGreaterEqual], (NormConeString)])0,{x,-1,1},{y,-1,1},{z,0,1}]

半正定値行列錐不等式

“SemidefiniteCone” ((SemidefiniteConeList[n])) 対称半正定値行列錐 x∈(ℝ^(n×n))

{{3,1},{1,4}}(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])IdentityMatrix[2]

N[Eigenvalues[{{3,1},{1,4}}-IdentityMatrix[2]]]

Quiet[RegionPlot[{{x,1},{1,y}}(Underscript[\[VectorGreaterEqual], (SemidefiniteConeString)])0,{x,0,2},{y,0,2}]]

RegionPlot3D[x z-y^2≥0,{x,0,4},{y,-3,3},{z,0,4},PlotPoints->50]

凸関数，凸集合

集合 C は，θ ∈ℝとCの中のx,y に対して θ x+(1-θ) y∈C であれば，アフィン(affine)集合とよばれる．

0≤ θ ≤ 1に対して上と同じことが成り立てば，Cは凸集合．

関数 f は，0≤ θ ≤ 1に対して

f (θx+(1-θ) y)≤f θ x+f (1-θ) y

が成り立てば，凸関数. (x, y ∈ dom (f) )

エピグラフは，f(x)≤ tとなる {x,t} の集合．

f[x_]:=x^2;
constrain[{x_,_}] := {x,f[x]};
plot = Plot[f[x],{x,-1,1},Filling->Top];
DynamicModule[{px={-1,1},py={1,1}},Show[plot,Graphics[{Dynamic[Line[{px,py}]],Locator[Dynamic[px,Function[px = constrain[#]]]],Locator[Dynamic[py,Function[py = constrain[#]]]]}]]]

凸最適化

凸最適化とは

凸集合上の凸関数の最小化．

凸最適化問題において，局所的な最小値がみつかれば，それが大域的な最小値であることが保証される．

凸最適化問題の標準形

 ({{minimize, f(x)}, {subject to, ({{g(x)≤0 }, {h(x)==0}})}}) 
 
ここで， f, g={(Subscript[g, 1]),…}  は凸関数，h はアフィン(affine) (i.e., h(x)== (Subscript[a, 0]).x+(Subscript[b, 0]))

“エピグラフ変換”を使うと, 上の問題は

 ({{minimize, t}, {subject to, ({{f(x)≤t }, {g(x)\[VectorLessEqual]0}, {h(x)==0}})}})

と書き直せる．つまり，どんな凸最適化問題も凸集合の上の線形関数の最適化問題として表現できる．

凸最適化の例

標準形

 ({{minimize, f(x)}, {subject to, ({{(Subscript[g, i])(x)≤0 }, {h(x)==0}})}}) 	f, (Subscript[g, i]): convex functions, h: affine

線形最適化 : 		(Subscript[f, 0])(x)=c.x,                                 (Subscript[g, i])(x)=(Subscript[a, i]).x+(Subscript[b, i])≥0

線形分数最適化：　(Subscript[f, 0])(x)=(α.x+β)/(γ.x+δ),  (Subscript[g, i])(x)=(Subscript[a, i]).x+(Subscript[b, i])≥0.

2次最適化：		(Subscript[f, 0])(x)=(1/2)x.Q.x+c.x,         (Subscript[g, i])(x)=(Subscript[a, i]).x+(Subscript[b, i])≥0.

2次錐最適化：		(Subscript[f, 0])(x)=c.x,                                 (Subscript[g, i])(x)=||(Subscript[A, i]).x+(Subscript[b, i])(Subscript[||, 2])≤(Subscript[α, i])x+(Subscript[β, i]).

モデル変換

最適化モデル(目的関数，制約条件，変数)を，ソルバが利用できるようにうまく変換する．

簡単な例で見てみよう．

||x(Subscript[||, 1])を最小化する問題の場合，新しいベクトル変数 s を(Subscript[s, i])==(Abs[(Subscript[x, i])])として導入し，制約条件

sconstraints = {s\[VectorGreaterEqual]0,-s\[VectorLessEqual]x\[VectorLessEqual]s}; (* ベクトル不等式! *)

を加えることで目的関数を線形関数に変換できる．

a=(({{-1, 1}, {1, 3}, {2, 1}})), b=(({{1}, {2}, {4}}))のもとで， ||a.x-b(Subscript[||, 1]) を最小化するx=(({{(Subscript[x, 1])}, {(Subscript[x, 2])}}))をみつける(2変数関数のフィッティングに相当)：

LinearOptimization[Total[s],{a=={{-1,1},{1,3},{2,1}},b=={1,2,4},-s\[VectorLessEqual]a.x-b\[VectorLessEqual]s},{x,s}]

||a.x-b((Subscript[||, 2])^2) +λ||x(Subscript[||, 1])を最小化するなら(QuadraticOptimizationを使って)

QuadraticOptimization[((Norm[a.x-b])^2)+λ Total[s], {a=={{-1,1},{1,3},{2,1}},b=={1,2,4},λ==3,sconstraints}, {x,s}]

となる．

モデル変換において重要な点は, 凸性を維持すること

よくある変換の例

・目的関数：||a.x-b(Subscript[||, 1])   ⇨   ({{minimize, ∑(Subscript[s, i])}, {subject to, ({{(Subscript[s, i])=|(Subscript[a, i]).x-(Subscript[b, i])| }, {-(Subscript[s, i])≤ (Subscript[a, i]).x-(Subscript[b, i])≤ (Subscript[s, i])}})}})

・目的関数：||a.x-b(Subscript[||, ∞])   ⇨    ({{minimize, s}, {subject to, ({{s=||a.x-b(Subscript[||, ∞])=(max_i)|(Subscript[a, i]).x-(Subscript[b, i])| }, {-(Subscript[s, i])≤ (Subscript[a, i]).x-(Subscript[b, i])≤ (Subscript[s, i])}})}})

・目的関数：|x|   ⇨   ({{minimize, (SuperPlus[x])+(SuperMinus[x])}, {subject to, ({{(SuperPlus[x])-(SuperMinus[x])=x}, {(SuperPlus[x])≥ 0, (SuperMinus[x])≥0}})}})

・目的関数：(max_i)((Subscript[a, i]).x-(Subscript[b, i]))   ⇨   ({{minimize, s}, {subject to, ({{s=(max_i)((Subscript[a, i]).x-(Subscript[b, i]))}, {(Subscript[a, i]).x-(Subscript[b, i])≤s}})}})

・目的関数：(min_i)((Subscript[a, i]).x-(Subscript[b, i]))　⇨　({{minimize, -s}, {subject to, ({{s=(min_i)((Subscript[a, i]).x-(Subscript[b, i]))}, {(Subscript[a, i]).x-(Subscript[b, i])≥s}})}})

線形最適化(Linear Optimization)

({{minimize, c.x}, {subject to, ({{a.x+b\[VectorGreaterEqual]0}, {(Subscript[a, eq]).x+(Subscript[b, eq])=0}})}})

・線形計画問題ともよばれる．

・機械精度，任意精度での計算を扱う専用ソルバ利用可能

LinearOptimization[objective, constraints, vars]
LinearOptimization[c,{a,b}]
LinearOptimization[c,{a,b},{aeq,beq}]

例

({{minimize, x+y}, {subject to, ({{x+2y-3≥ 0}, {x≥-1, y≥-1}})}})

LinearOptimization[x+y,{x+2y≥3,x≥-1,y≥-1},{x,y}]

LinearOptimization[{1,1},{{{1,2},{1,0},{0,1}},{-3,1,1}}]

以下，いくつかの等価な入力方法

ベクトル変数v={x,y}を用いて

LinearOptimization[{1,1}.v,{{1,2}.v≥3,{1,0}.v≥-1,{0,1}.v≥-1},v]

VectorGreaterEqual (or VectorLessEqual)を用いて複数の制約条件をまとめる:

LinearOptimization[{1,1}.v,VectorGreaterEqual[{{{1,2},{1,0},{0,1}}.v,{3,-1,-1}}],v]

Esc v>= Esc で VectorGreaterEqual をコンパクトに “ \[VectorGreaterEqual] ” と表記:

LinearOptimization[{1,1}.v,{{1,2},{1,0},{0,1}}.v\[VectorGreaterEqual]{3,-1,-1},v]

パラメータを指定する等式を利用:

LinearOptimization[c.v,{a.v\[VectorGreaterEqual]b,c=={1,1},a=={{1,2},{1,0},{0,1}},b=={3,-1,-1}},v]

線形最適化: ロバストなデータフィッティング

データ: {{(Subscript[x, 1]),(Subscript[y, 1])},{(Subscript[x, 2]),(Subscript[y, 2])},...,{(Subscript[x, n]),(Subscript[y, n])}}  に対し，||a.λ-b(Subscript[||, 1])を最小化することで，ロバストな非線形フィットを行う

a=.;b=.

data=Block[{d,xOut},d=Table[{x,Sin[2 x]+Cos[(x^2)]+RandomVariate[NormalDistribution[0,0.5]]},{x,-3,3,0.04}];xOut=RandomSample[Range[-3,3,0.1],30];Join[d,Transpose[{xOut,(1/2) (Sin[2 xOut]+Cos[(xOut^2)])+RandomReal[{3,6},Length[xOut]]}]]];
dataPlot = ListPlot[data, PlotRange->All]

基底関数として (Subscript[ϕ, i])(x)=(((1+(((x-(Subscript[c, i])))^2)))^(1/2)) を用いる. フィッティング関数はその線形結合 s(x)=(∑_(i=1))(Subscript[λ, i])(Subscript[ϕ, i])(x). 
データから計画行列を構成．

bases = Sqrt[1+(x-#)^2]&/@Subdivide[-3,3,10];
input=DesignMatrix[data,bases,x,IncludeConstantBasis->False];

Plot[bases,{x,-3,3}]

各基底の係数 (Subscript[λ, i]) を決めたい. 
モデル変換 ||a.λ-b(Subscript[||, 1])=(∑_i)(Subscript[s, i]),(Subscript[s, i])=|(Subscript[a, i]).λ-(Subscript[b, i])|, -(Subscript[s, i])≤ (Subscript[a, i]).λ-(Subscript[b, i])≤ (Subscript[s, i]) (i=1,2,..,n) を使う
(a: input, b: output, sは各成分の偏差)

res=LinearOptimization[Total[s],{a==input,b==data[[All,2]],-s\[VectorLessEqual]a.λ-b\[VectorLessEqual]s},{λ,s}];
Short[res,6]

こうして見つかったフィッティングカーブは外れ値に対して強い(ロバスト)ことがわかる:

Show[dataPlot,Plot[bases.λ/.res,{x,-3,3}]]

最小二乗法との比較

lsqfit=Fit[data,bases,x];

Plot[{Sin[2x]+Cos[x^2],Evaluate[bases.λ/.res],lsqfit},{x,-3,3}, PlotLegends->{"Reference","(Subscript[L, 1])フィット","最小二乗法"}]

線形分数最適化(Linear Fractional Optimization)

({{minimize, (α.x+β)/(γ.x+δ) }, {subject to, ({{a.x+b\[VectorGreaterEqual]0}, {(Subscript[a, eq]).x+(Subscript[b, eq])=0}})}})

LinearFractionalOptimization[objective, constraints, vars]
LinearFractionalOptimization[{α,β,γ,δ},{a,b}]
LinearFractionalOptimization[{α,β,γ,δ},{a,b},{aeq,beq}]

・γ = 0 なら線形最適化と同じ.
・y = (x/(γ x + δ)), z=(1/(γ x + δ))=>x=y/zと変換することで，線形最適化問題

({{minimize, α.y+β t}, {subject to, ({{a.y+b z\[VectorGreaterEqual]0}, {(Subscript[a, eq]).y+(Subscript[b, eq])z=0}})}})

として計算可能．

・コストを最小化，利益を最大化，などの複数の目的関数の同時最適化問題に適用できる．

線形分数最適化: 輸送問題

5種類の品物を輸送する際の輸送コストを最小限にしつつ，利益を最大化したい．

品物のそれぞれを1単位輸送するのにかかるコストと利益

costPerUnit = {100,100,100,100,100};
profitPerUnit={115,108,130,120,110};

それぞれの重量と，積載可能な単位数

weightPerUnit={100,100,110,200,115};
maxUnits={30,20,20,30,20};

トラックを利用する際の最小コストは，1回あたり$2000(人件費等).  全輸送コストは
(xは各品の単位数を成分とするベクトル)

cost=2000+costPerUnit.x;

トラック1回の輸送で，輸送量によらず$1000の利益があるとする．全利益は

profit=1000+profitPerUnit.x;

トラック1台あたりの最大積載量は8000ポンド．積載量に関する制約条件は

wtConstraint=weightPerUnit.x≤ 8000;
goodsConstraints=0\[VectorLessEqual]x\[VectorLessEqual]maxUnits;

最適な品物数の組み合わせは，輸送コストと利益の比(cost/profit)を最小化することで得られる．
(cost/profitはefficiencyともよばれる)

res = LinearFractionalOptimization[cost/profit,{wtConstraint,goodsConstraints},x∈Vectors[5]]

TableForm[{x/.res}, TableHeadings->{{"Units"},Array[(Subscript[x, #])&,5]}]

二次最適化(Quadratic Optimization)

({{minimize, (1/2)xᵀ.Q.x + c.x }, {subject to, ({{a.x+b\[VectorGreaterEqual]0}, {(Subscript[a, eq]).x+(Subscript[b, eq])=0}})}})

QuadraticOptimization[objective, constraints, vars]
QuadraticOptimization[{Q,c},{a,b}]
QuadraticOptimization[{Q,c},{a,b},{aeq,beq}]]
QuadraticOptimization[{{Qf,0},c},...]
QuadraticOptimization[{{Qf,cf},c},...]

・行列Qは正値対称行列．

・Q=0なら線形最適化問題と等価．

・2次形式の部分は，積に分けた形，xᵀ.(Subscript[Q, f])ᵀ(Subscript[Q, f]).x，あるいは (xᵀ.(Subscript[Q, f])ᵀ+(Subscript[c, f])).((Subscript[Q, f]).x+(Subscript[c, f])) でもOK．

例1

({{minimize, 2(x^2)+20(y^2)+6x y+5x }, {subject to, -x+y-2≥ 0}})

obj=2(x^2)+20(y^2)+6x y+5x;
res=QuadraticOptimization[obj,-x+y≥2,{x,y}]

Show[Plot3D[obj,{x,-2,0},{y,0,1},RegionFunction->Function[{x,y,z},-x+y≥2],MeshFunctions->{#3&},Mesh->20,BoxRatios->1,ImageSize->300,AxesLabel->{Style[x,Large],Style[y,Large],Style["obj",Large]}],Graphics3D[{Red, PointSize[0.05],Point[{x,y,obj}/.res]}]]

例2

({{minimize, (x^2)+(y^2)}, {subject to, ({{0≤y≤2}, {x+2y=4}})}})

res=QuadraticOptimization[(x^2)+(y^2),{x+2y==4,0≤y≤2},{x,y}]

目的関数を1/2x.q.x+c.x として，制約条件を a.x+b\[VectorGreaterEqual]0および (Subscript[a, eq]).x+(Subscript[b, eq])==0として定義する：

{q,c}={{{2,0},{0,2}},{0,0}};
{a,b}={{{0,1},{0,-1}},{0,2}};
{aeq,beq}= {{{1,2}},{-4}};

行列ベクトル入力を使って解く：

QuadraticOptimization[{q,c},{a,b},{aeq,beq}]

最適点は (x^2)+(y^2)の等高線が等式制約に接するところにある：

Show[RegionPlot[Evaluate[(x^2)+(y^2)≤((x^2)+(y^2)/.res)],{x,0,2},{y,0,2},Epilog->{Blue,PointSize[0.04],Point[{x,y}/. res]},ImageSize->200],ContourPlot[x+2y==4,{x,0,2},{y,0,2}, ContourStyle->Black]]

例3

a=.;b=.;

3 次元データの集合P, Q (set1, set2)を分離する2次曲面をみつける

{set1,set2}=Block[{},{RandomReal[{-1,1},{200,3}]/. {x_,y_,z_}:>{x,-2 (x^2)+2 (y^2),z},Transpose[{0,-1,0}+Transpose[RandomReal[{-(1/2),(1/2)},{200,3}]]]}];
dataplot=ListPointPlot3D[{set1,set2}, PlotRange->All,PlotStyle->{{Red,PointSize[0.025]},{Blue,PointSize[0.025]}},AspectRatio->1]

データ点情報を基底関数を通して並べた計画行列を構築

{data1,data2}=DesignMatrix[#,{x,(x^2),y,(y^2),x y},{x,y}]&/@{set1,set2};

set1, set2を分離するにはset1の点はa.(Subscript[p, i])+b≥1, set2の点はa.(Subscript[p, i])+b≤ -1 を満たす必要がある．(制約条件)
最適な分離超曲面を||a((||)^2)を最小化することで見つける

res=QuadraticOptimization[((Norm[a])^2),{data1.a+b\[VectorGreaterEqual]1,data2.a+b\[VectorLessEqual]-1},{a,b},Method->"COIN"]

分離超曲面を表す多項式は

poly=(a.{1,x,(x^2),y,(y^2),x y}+b)/.res

データ点とともに描画してみる

Show[dataplot, ContourPlot3D[poly==0,{x,-1.5,1.5},{y,-1.5,1.5},{z,-1,1},Mesh->None,ContourStyle->Directive[Orange,Opacity[0.6],Specularity[White,30]],Lighting->"Neutral"]]

錐最適化(Conic Optimization)

錐 (cone)

ベクトル空間Vの部分空間Cをとり, Cの要素xをα倍 (α > 0) したαxもCの要素であるようなとき，Cは錐であるという．

凸最適化の標準形：  ({{minimize, c.x}, {subject to, ({{h(x)==0}, {(Subscript[g, j])(x)≤0}})}}), j=1,…,k

において，h(x)がアフィン変換，(Subscript[g, j])(x)≤0 (j=1,…,k)が錐(cone)を表す場合，

({{minimize, c.x}, {subject to, ({{(Subscript[a, 0]).x+(Subscript[b, 0])==0}, {(Subscript[a, j]).x+(Subscript[b, j])(Underscript[\[VectorGreaterEqual], (Subscript[κ, j])])0}})}})

錐線形最適化問題とよばれる．ここで，(Subscript[κ, j])は錐体(凸集合)．

凸関数 (Subscript[f, i]) に対して，条件 (Subscript[f, i])(x) ≤ 0 が (Subscript[a, i]).x+(Subscript[b, i])(Underscript[\[VectorGreaterEqual], (Subscript[κ, i])])0 と等価となるような凸錐 (Subscript[κ, i]) が存在する[http://www.seas.ucla.edu/~vandenbe/236C/lectures/conic.pdf] .

Mathematicaの最適化ソルバでは以下の錐体が指定可能:

これだけあれば，多くの問題を扱うのに十分．

モデル変換 revisited

凸モデル変換の一般的な目的は，与えられた問題を錐線形最適化問題に変換すること．

たとえば，フィッティング問題
	((Norm[a.x-b])^2)+λ Norm[x,1] を最小化
は，錐線形最適化問題に帰着する．第1項，第2項ともに凸なので，

({{minimize, s + λ t}, {subject to, ({{(u^2)=s}, {Norm[a.x+b]≤u}, {Norm[x,1]≤t}})}})

とできて，これは以下のように解ける:

parameters = {a=={{1,2},{3,4},{5,6}},b=={7,8,9},λ==1};

ConicOptimization[s+λ t,{u^2≤s,Norm[a.x-b]≤u,Norm[x,1]≤t,parameters},{x,s,t,u}]

内部では，ConicOptimization関数が錐線形問題の標準形に変換している:

conicConstraints = ConicOptimization[s+3t,{u^2≤s,Norm[a.x+b]≤u,Norm[x,1]≤t,parameters},{x,s,t,u},"ConicConstraints"]

conicConstraints[[All,1,1,1,1]]=Map[MatrixForm,conicConstraints[[All,1,1,1,1]]];
conicConstraints

感度解析

利益を最大化する製品(3種類)の組み合わせを見つけたい． 各製品は4つの機械を使って製造する．それぞれの製品製造に各機械が必要とする時間は以下の通り．また，各機械は1週間に最大2400分稼働可能とする．

machineTime ={{20,10,10},{12,28,16},{15,6,5},{10,15,0}};Row[BarChart[machineTime[[#]],LabelingFunction->"RadialCenter",PlotLabel->"machine StandardForm[#1]",ChartLabels->{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])},ChartElementFunction->"GlassRectangle",ChartStyle->101,BarSpacing->None,ImageSize->200]&/@Range[4]]

製品(Subscript[x, 1]), (Subscript[x, 2]), (Subscript[x, 3])は最大製造量はそれぞれ，100単位, 40単位, 60単位である．
それぞれ1単位で得られる利益は，$45, $60, $50．

timeConstraint=machineTime.x\[VectorLessEqual]2400;  (* 4-dim *)
capacity=x\[VectorLessEqual]{100,40,60}; (* 3-dim *)
profit = {45.,60.,50.}.x;

以上の条件下での最大の利益は(“PrimalMinimumValue”指定で目的関数の最適値出力)

max0 = -LinearOptimization[-profit,{timeConstraint,capacity},x,"PrimalMinimumValue"]

それを実現する製造量は

xmin= LinearOptimization[-profit,{timeConstraint,capacity},x,"PrimalMinimizerVector"]

感度解析

最大製造量を変化させたときの利益に対する影響を見たい
“ConstraintSensitivity”を指定するとbに対する目的関数の変化率{∂(SuperStar[f])/∂b,∂(SuperStar[f])/∂(Subscript[b, eq])}を返す．

{sens,{a,b}}=LinearOptimization[-profit,{timeConstraint,capacity} ,x,
{"ConstraintSensitivity","LinearInequalityConstraints"}]

(最後3つの数0, -7.1(x10^(-15)), -10.5がcapacityに対する感度)

b, (Subscript[b, eq])の変化に対して感度がある条件をピックアップ

spos=Flatten@Position[Chop[sens[[1]]],_?(Abs[#]>0&)];
Table[a[[i]].{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}+b[[i]]≥0,{i,spos}]//TableForm

製品3の最大製造可能量を増やすと，利益が増す．

x1, x2はサチってないため現制約条件の値から変えても変化ない．

max1=-LinearOptimization[-profit,{timeConstraint,x\[VectorLessEqual]{100,40,60+10}} ,x,"PrimalMinimumValue"]

この増分は，上の感度解析から得られる値と一致する．

max1-max0

-sens[[1,spos[[3]]]]×10

Summary

凸最適化は，Wolfram言語環境の中で高度に統合された強力なツール．

求解可能な問題のクラスが大幅に拡大

最適解の詳細情報を即座に提供するソルバ

Wolfram言語に新たに加わったモデリングツール

一般化されたベクトル不等式

(自動)次元推定

パラメータの特定法(制約条件)