第1回輪講 2011.5.10 (火)

慶應義塾大学理工学部
管理工学科４年曹研究室
60803571
遠藤健司
 凸集合、凸関数
 非線形計画問題を解く
 プログラミング

「数理計画法の基礎」



「数理計画法」～最適化の手法



坂和正敏著
森北出版（株）（1999,初版）
森哲夫著
共立出版（株）（1994,初版
2010年度「管理工学用数学第２」授業用テキスト
def
1
2
n
1
2
0    1}
X  R が凸集合  [x , x ]  {x  R | x  x  (1   )x , n
: X に属する任意の２点 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 )
“＞”の場合：強意凸関数
“≦”の場合：凹関数
定理：１回微分可能な凸関数
１回連続微分可能な関数 f (x) が凸関数であるための必要十分条件は、
1
2
n
1
2
2
1
2
任意の x , x  R に対して、 f (x )  f (x )  f (x )(x  x ) が成り立つこ
とである。
定理：２回微分可能な凸関数
２回連続微分可能な関数 f (x) が凸関数であるための必要十分条件は、
f (x) のヘッセ行列 2 f (x) が任意の x  R n に対して半正定である
ことである。（ yT 2 f (x)y  0, y  Rn ）
また、２回連続微分可能な関数 f (x) が強意凸関数であるための十分
n
条件は、 f (x) のヘッセ行列 2 f (x) が任意の x  R に対して正定
であることである。（ yT 2 f (x)y  0, y  Rn , y  0 ）
ヘッセ行列が正定
⇔ ヘッセ行列の任意
の首座小行列式＞０
 2
f ( x)

2

x
 1

2
 f (x) :  ・・・

 2
f ( x)

 xn x1
・・・
2
f ( x)
xi x j
・・・

2
f ( x) 
x1xn


・・・



2
f ( x) 
2
xn


3変数2次関数の非線形計画問題（NLP）
min f (x)  3x1  x2  x3  2 x1 x2  x2 x3  3x1 x3  4 x1  6 x2  3x3  6
2
2
2
s.t.　g (x)  x2  3x3  5  0
2

ラグランジュ緩和問題
min L(x,  )  3 x1  x2  x3  2 x1 x2  x2 x3  3 x1 x3  4 x1  6 x2  3 x3  6
2
2
2
  ( x2  3 x3  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)
x1x2
2
f ( x)
2
x2
2
f ( x)
x3x2
2  3
 6


  2
2 1
  3 1 2 



2
f ( x) 
x1x3

2


f ( x) 
x2 x3

2

f ( x) 
2
x3

2
ヘッセ行列  f (x) の首座小行列について
6
2
3
6 2
| 6 | 6  0,  8  0, 2
2 1  4  0
2 2
 3 1 2
ヘッセ行列  f (x) の任意の首座小行列が０よりも大きいため、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
１．初期点 xl を選び、l : 1 とする。
２．現在の点 x l において停止基準をみたせば終了、そうでない場合は降下方向
dl  T f (xl ) を求める。
f (xl )dl
l
, l  1,2,...
３．一次元探索問題を解き、ステップ幅    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

e  g
  
2b f    h 
f 2c   i 
d








































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]+"))");
降下方向ベクトル dl  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
と分母部分 (dl )T Qdl を求め、最適ス
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 dl を求める。
f (xl 1 ) を求める。
ここで、探索過程で導出され
た f を配列に格納していく
停止基準より、|| f (x )  f (x ) || 
ならば探索を終了する。
今回、  1.0e7 とした。
l 1
l
定理：凸関数の最適性の必要十分条件
T
2
0
n
d
n
  f (x )d  0, d  R
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  3x3  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  3x1 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 gi (x0 )  0
0
0
0
0
i 1
gi (x 0 )  0, i gi (x 0 )  0, i  0, i  1,2,...,m
0
0
 x L(x ,  )  f (x )  g (x )
 6 x1  2 x2  3x3  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

Download Report