非線形計画問題 自作 - 慶應義塾大学 理工学部管理工学科
Download
Report
Transcript 非線形計画問題 自作 - 慶應義塾大学 理工学部管理工学科
慶應義塾大学 理工学部
管理工学科4年 曹研究室
60803571
遠藤 健司
凸集合、凸関数
非線形計画問題を解く
プログラミング
「数理計画法の基礎」
「数理計画法」~最適化の手法
坂和正敏 著
森北出版(株)(1999,初版)
森哲夫 著
共立出版(株)(1994,初版
2010年度「管理工学用数学第2」授業用テキスト
def
X R が凸集合
n
[x1 , x 2 ] {x R n | x x1 (1 )x 2 , 0 1}
: X に属する任意の2点 x1 , x 2 を結ぶ線分が必ず X に含まれるとき、
集合 X は凸集合である。
x1
x1
x2
x1
x1
x2
x2
x2
凸集合
非凸集合
関数のグラフの「上側」の部分(エピグラフ)が凸集合になっている関数
=下に膨らんだ形をする関数
def
f が凸関数
x1 , x 2 R n , 0 1
f (x1 ) (1 ) f (x 2 ) f (x1 (1 )x 2 )
“>”の場合:強意凸関数
“≦”の場合:凹関数
定理:1回微分可能な凸関数
1回連続微分可能な関数 f (x) が凸関数であるための必要十分条件は、
1
2
n
1
2
2
1
2
任意の x , x R に対して、 f (x ) f (x ) f (x )(x x ) が成り立つこ
とである。
定理:2回微分可能な凸関数
2回連続微分可能な関数 f (x) が凸関数であるための必要十分条件は、
f (x) のヘッセ行列 2 f (x) が任意の x R n に対して半正定である
ことである。( yT 2 f (x)y 0, y R n )
また、2回連続微分可能な関数 f (x) が強意凸関数であるための十分
n
条件は、 f (x) のヘッセ行列 2 f (x) が任意の x R に対して正定
であることである。( yT 2 f (x)y 0, y R n , y 0 )
ヘッセ行列が正定
⇔ ヘッセ行列の任意
の首座小行列式>0
2
f ( x)
2
x
1
2
f (x) : ・・・
2
f ( x)
x
x
n
1
・・・
2
f ( x)
xi x j
・・・
2
f ( x)
x1xn
・・・
2
f ( x)
2
xn
3変数2次関数の非線形計画問題(NLP)
min
f (x) 3x1 x2 x3 2 x1 x2 x2 x3 3 x1 x3 4 x1 6 x2 3x3 6
2
2
2
s.t. g (x) x2 3x3 5 0
2
min
ラグランジュ緩和問題
L(x, ) 3x1 x2 x3 2 x1 x2 x2 x3 3x1 x3 4 x1 6 x2 3x3 6
2
2
2
( x2 3x3 5)
2
s.t.
0
f ( x ) (
f (x) f (x) f (x) T
,
,
)
x1
x2
x3
6 x1 2 x2 3 x3 4
2 x1 2 x2 x3 6
3x x 2 x 3
1
2
3
2
f ( x)
2
x
1
2
2
f ( x)
f ( x)
x
x
22 1
x x f (x)
3 1
2
f ( x)
x1x2
2
f ( x)
2
x2
2
f ( x)
x3x2
2 3
6
2
2 1
3 1 2
2
f ( x)
x1x3
2
f ( x)
x2 x3
2
f ( x)
2
x3
2
ヘッセ行列 f (x) の首座小行列について
| 6 | 6 0,
6 2
2 2
8 0,
6
2
3
2
2
1 4 0
3 1
2
ヘッセ行列 f (x) の任意の首座小行列が0よりも大きいため、f (x)
は正定となる。
よって、 f (x) は強意凸関数であることがわかった。
また、 g (x ) は明らかに下に凸な形のグラフであるため、 g (x ) は凸関
数であることがわかる。
2
では、ラグランジュ緩和問題の局所的最適解を、最急降下法を用いて
求める。
→JAVAを用いた最急降下法のプラグラムの作成
自力でプログラムを作成した結果、緩和問題の最急降下法プロ
グラミングはできなかった…。
しかし、制約条件なしの非線形最適化問題(NLP)を最急降下法で解くこ
とのできるプログラミングの作成に成功!
ということで、話が少々逸脱しますが、JAVAで作成したNLPの最急降下法
プログラミングをお見せします。
NLP6.java
min
f (x) ax1 bx2 cx3 dx1 x2 ex2 x3 fx1 x3 gx1 hx2 ix3 j
2
2
2
3x1 x2 x3 2 x1 x2 x2 x3 3x1 x3 4 x1 6 x2 3x3 6
2
2
2
→ Wolfram Mathematica 7.0
1.初期点 xl を選び、l : 1 とする。
2.現在の点 x l において停止基準をみたせば終了、そうでない場合は降下方向
d l T f (xl ) を求める。
f (xl )d l
l
,
l 1,2,...
3.一次元探索問題を解き、ステップ幅 l T
(d ) Qd l
l 1
l
l l
を求め、 x : x d , l : l 1
として、手順2.へ戻る。
停止基準について
l
1. || f (x ) || となれば終了。
2.関数の変化量がある許容値以内
ならば終了。
|| f (xl 1 ) f (xl ) ||
f (x) ax1 bx2 cx3 dx1 x2 ex2 x3 fx1 x3 gx1 hx2 ix3 j
2
2
2
T
f g
2a d
1 T
1 T
T
x Qx b x j x d 2b e x h x j
2
2
e 2c i
f
f ( x ) x T Q b
2a
xT d
e
d
2b
f
e g
f h
2c i
import java.io.*;
public class NLP6
{
public static void main (String[] args) throws IOException
{
double a,b,c,d,e,f,g,h,i,j,x1s,x2s,x3s,QQ,DD,A,AU,AD,F,FF,fx;
int m,n;
int p=1;
int t=0;
String s;
double x[] = new double[3];
double Q[][] = new double[3][3];
double B[] = new double[3];
double Q0[] = new double[3];
double D[] = new double[3];
double D0[] = new double[3];
double xs[] = new double[3];
double F0[] = new double[3];
double z[] = new double[100000];
InputStreamReader in = new InputStreamReader(System.in);
BufferedReader br = new BufferedReader(in);
System.out.println("min f(x) =
a*x1^2+b*x2^2+c*x3^2+d*x1x2+e*x2x3+f*x3x1+g*x1+h*x2+i*x3+j");
System.out.print("a=");
s = br.readLine();
a = Double.valueOf(s).doubleValue();
System.out.print("b=");
s = br.readLine();
b = Double.valueOf(s).doubleValue();
System.out.print("c=");
s = br.readLine();
c = Double.valueOf(s).doubleValue();
System.out.print("d=");
s = br.readLine();
d = Double.valueOf(s).doubleValue();
System.out.print("e=");
s = br.readLine();
e = Double.valueOf(s).doubleValue();
System.out.print("f=");
s = br.readLine();
f = Double.valueOf(s).doubleValue();
変数と配列(ベクトル)の定義
非線形計画問題の係数(a~j)
の入力値を取り込む
System.out.print("g=");
s = br.readLine();
g = Double.valueOf(s).doubleValue();
System.out.print("h=");
s = br.readLine();
h = Double.valueOf(s).doubleValue();
System.out.print("i=");
s = br.readLine();
i = Double.valueOf(s).doubleValue();
System.out.print("j=");
s = br.readLine();
j = Double.valueOf(s).doubleValue();
System.out.println("初期点(x1s,x2s,x3s)を入力してください");
System.out.print("x1s=");
s = br.readLine();
x1s = Double.valueOf(s).doubleValue();
System.out.print("x2s=");
s = br.readLine();
x2s = Double.valueOf(s).doubleValue();
System.out.print("x3s=");
s = br.readLine();
x3s = Double.valueOf(s).doubleValue();
x[0]=x1s;
x[1]=x2s;
x[2]=x3s;
Q[0][0]=2*a;
Q[0][1]=d;
Q[0][2]=f;
Q[1][0]=d;
Q[1][1]=2*b;
Q[1][2]=e;
Q[2][0]=f;
Q[2][1]=e;
Q[2][2]=2*c;
B[0]=g;
B[1]=h;
B[2]=i;
NLPの係数(a~j)の入力値を
取り込む
初期値 xl の入力値を取り込む
取り込んだNLPの係数と変数の初期値の各々を
配列へと格納していく
do
{
for(m=0; m<3; m++)
{
QQ=0;
for(n=0;n<3;n++)
{
QQ = QQ + Q[m][n]*x[n];
Q0[m] = QQ + B[m];
}
}
for(m=0; m<3; m++)
{
D[m] = Q0[m];
}
if(Math.abs(D[0])<1.0e-10 && Math.abs(D[1])<1.0e-10 &&
Math.abs(D[2])<1.0e-10)
{
break;
}
System.out.println("探索"+p+"回目の降下方向ベクトル
は,(d1,d2,d3)=(-("+D[0]+"),-("+D[1]+"),-("+D[2]+"))");
降下方向ベクトル d l T f (xl )
を求める
停止基準より、|| f (xl ) ||
ならば、探索を終了する
今回、 1.0e 10 とした。
AU=Q0[0]*D[0]+Q0[1]*D[1]+Q0[2]*D[2];
for(m=0; m<3; m++)
{
DD=0;
for(n=0;n<3;n++)
{
DD = DD + Q[m][n]*D[n];
D0[m] = DD;
}
}
AD=D[0]*D0[0]+D[1]*D0[1]+D[2]*D0[2];
A=AU/AD;
System.out.println("探索"+p+"回目の最適ステップ幅
は,a="+A);
最適ステップ幅の分子部分 f (x )d
と分母部分 (d l )T Qd l を求め、最適ス
l
テップ幅 を求める。
l
l
for(m=0; m<3; m++)
{
xs[m] = x[m]-A*D[m];
x[m] = xs[m];
}
System.out.println("探索"+p+"回目における近似解
は,(x1,x2,x3)=("+xs[0]+","+xs[1]+","+xs[2]+")");
for(m=0; m<3; m++)
{
FF=0;
for(n=0;n<3;n++)
{
FF = FF + Q[m][n]*x[n];
F0[m] = FF;
}
}
fx =
j+(B[0]*x[0]+B[1]*x[1]+B[2]*x[2])+0.5*(F0[0]*x[0]+F0[1]*x[1]+F0
[2]*x[2]);
t+=1;
z[t]=fx;
System.out.println("その時,fの値は, f(x)="+z[t]);
p+=1;
}while( Math.abs(z[t]-z[t-1])>1.0e-7 );
System.out.println(" ");
System.out.println("よって最適解
は,(x1*,x2*,x3*)=("+xs[0]+","+xs[1]+","+xs[2]+")");
System.out.println("最適値は, f(x*)="+z[t]);
}
}
ステップ幅によって導かれた点
xl 1 : xl l d l を求める。
f (x l 1 ) を求める。
ここで、探索過程で導出され
た f を配列に格納していく
l 1
l
停止基準より、|| f (x ) f (x ) ||
ならば探索を終了する。
今回、 1.0e 7 とした。
定理:凸関数の最適性の必要十分条件
T
2
0
n
d
d
R
n
f ( x )d 0,
f (x) が R 上で微分可能な凸関数とするとき、x が大域的最適解であるため
の必要十分条件は、 f (x ) 0 が成立することである。
3 11
x ( x1 , x2 , x3 ) ( , ,2)
2
2
3
11
6 * 2 * 3 * (2) 4
2
2
0
6 x1 2 x2 3 x3 4
3
11
f (x) 2 x1 2 x2 x3 6 2 * 2 * 1* (2) 6 0 0
2
2
0
3
x
x
2
x
3
1
2
3
3
11
3 * 2 1* 2 2 * (2) 3
f (x) は先述したように、強意凸関数であるので、x は大域的最適解となり、
最小値は f (x ) 16.5
では本題ということで、先程の3変数2次関数NLPを解く。
min
f (x) 3x1 x2 x3 2 x1 x2 x2 x3 3 x1 x3 4 x1 6 x2 3x3 6
2
2
2
s.t. g (x) x2 3x3 5 0
2
Wolfram Mathematica 7.0
FindMinimum関数を用いて解くと
7
4
x ( x1 , x2 , x3 ) ( ,3, )
3
3
86
f (x )
9
7
4
x ( x1 , x2 , x3 ) ( ,3, )
3
3
86
f (x )
9
定理:凸計画問題に対する最適性の十分条件
不等式制約のある非線形計画問題において、f (x), gi (x), i 1,..., m が全て凸関数
のとき(つまり凸計画問題のとき)、x X において、KKT条件を満たして
いれば、x は大域的最適解となる。
Kuhn-Tucker, KKT 条件
m
x L(x , λ ) f (x ) i g i (x 0 ) 0
0
0
0
0
i 1
g i (x 0 ) 0,
i 0 g i (x 0 ) 0,
i 0 0,
i 1,2,..., m
x L( x , ) f (x ) g (x )
6 x1 2 x2 3 x3 4
0
2 x1 2 x2 x3 6 2 x2
3x x 2 x 3
3
1
2
3
4
7
6 * 2 * ( 3) 3 * 4
3
0
0
0
3
7
4
10
2 * 2 * ( 3) * 6 2 * ( 3) 6 0
3
3
3
3
3
7
4
5
3 * *(3) 2 * 3
3
3
3
5
9
*
4
2
g (x ) x22 3x3 5 3 3 * 5 0
3
5
* g (x ) * 0 0
9
5
* 0
9
5
よって、KKT条件を満たすラグランジュ乗数 が存在し、KKT条件が
9
成立する。
7
4
今回のNLPは、先述の通り f (x), g (x) は凸関数であるため、x ( x , x , x ) ( 3 ,3, 3 )
86
が大域的最適解となり、最小値は f (x ) 9
*
1
2
3