非線形計画法

Ⅵ．非線形計画法

非線形計画法

　　今まで扱ってきたのは全て離散的なデータでした．離散的な世界の中で，データを記憶し，探索するためのデータ構造やアルゴリズムについて説明してきました．ここでは，連続的なデータによって表される世界について考えてみます．一般に，連続的なデータは関数によって表すことができます．従って，関数の形さえ記憶しておけば，個々のデータを記憶する必要はありませんので，記憶に要求される領域は非常に小さくなりますが，関数内のある点を探索しようとした場合，今まで述べてきたようなアルゴリズムをそのまま使用するわけにはいきません．

　　探索を目的とする点としては，関数の値を 0 とする点（方程式の根など），関数の値が最大（関数の符号を変えれば，最小）となる点，などが存在しますが，ここでは，関数（目的関数）の値が最大となる点を探索すること，つまり，最適化問題について考えます．一般に，関数に現れる変数の値は -∞ ～ +∞ となりますので，何らかの制約条件（これも，関数で表現される）が付加されます．ここでは，目的関数，または（及び），制約条件に非線形関数を含む場合，つまり，非線形計画法（ NP: Nonlinear Programming ）について説明します．非線形計画法が対象とする問題は以下のように記述できます．

目的関数，

y = f(x)　　x = [x₁, x₂, ･･･, x_n]^T　　　(1)

を，制約条件，

h_i(x) ＝ 0,　　i＝1，2，･･･，m　　　(2)

g_j(x) ≧ 0,　　j＝1，2，･･･，l　　　(3)

のもとで，最大（最小）にする．

　　非線形計画法において使用されるアルゴリズムは，目的関数の勾配情報（微分）を利用するか否かによって，大きく 2 つに分類できます．ここでは，勾配を利用する方法に対する基本アルゴリズムについて説明します．その基本アルゴリズムは以下のようになります．
1. 適当な x の初期値 x(0) を与え，繰り返し回数を表すパラメータ k を 0 にしておきます．
2. x(k) が解であるか否かの判定を行います．つまり，最小点であるための必要条件（勾配ベクトル = 0 ），
  
  　　　∥∇f(x(k))∥ = ∥[∂f/∂x₁，･･･，∂f/∂x_n]^T∥ = 0
  
  を調べることによって行います．実際上は，厳密に 0 になることはほとんどありませんので，上記の値が，適当な小さな値（ε₁）より小さいか否かによって行います．
  
  　　上記の判定だけですと，場合によっては，プログラムが無限ループに入ってしまう可能性もあります．従って，繰り返し回数がある回数以上になったら終了させることはもちろん，以下に示すような判定基準を併用することもあります．
  
  ∥f(x(k+1)) - f(x(k))∥ < ε₂
  
  ∥x(k+1) - x(k)∥ < ε₃
  
  ∥f(x(k))∥ < ε₄
3. 上記の判定の結果，解に収束していない場合は，k 回目の計算結果を基にして，関数値 f を減少させる新しい点 x(k+1) を生成し（右図参照），k = k + 1 として，ステップ ⅱ へ戻ります．新しい点 x(k+1) は，式，
  
  　　　x(k+1) = x(k) + α(k)p(k)
  
  によって計算されます．そのためには，探索方向 p(k) とその方向に沿った歩み幅（ step size ） α(k) を決める必要があります．
  
  　　探索方向 p(k) を決めるための情報として，勾配ベクトル ∇f(x(k)) = [∂f/∂x₁，･･･，∂f/∂x_n]^T や Hesse 行列
  
  等が利用されます．この p(k) の計算方法によって，アルゴリズムが異なってきます．一般に，多くの情報を用いれば繰り返しの回数は少なくなりますが，その反面，1 回の点を計算するために必要な計算量が増加します．
  
  　　　歩み幅 α(k) を一定値に固定すると，一般に，f の減少が保証できなくなり，その結果アルゴリズムの収束性が成立しなくなる場合があります．そこで，1 次元探索によって，つまり，
  
  f(x(k) + α(k)p(k))
  
  を最小にする α(k) を選択する方法がよく使用されます．
　　以上述べたアルゴリズムにおいて，最も重要な問題は，局所的最適解の問題です．先の図に示したように，山（今まで，最小値を求める方法について説明しましたので，「くぼみ」といった方が良いかもしれませんが，ここでは，最大値を求める問題にすり替え，山という表現を使用します）が一つだけの場合は問題ありませんが，一般には，右図に示すように，多くの山が存在します（多峰性）．最適化アルゴリズムは，これらの山の中から，最も最適なもの（最も高い山）を選び出す必要があります．しかし，上記のアルゴリズムを適用した場合，収束先は，ほとんど初期値によって決まってしまいます．例えば，右図の場合，点 A を初期値とすれば左の山（最適解）へ，また，点 B を初期値とすれば右の山（局所的最適解）へ収束してしまう可能性が非常に高いことになります．

最急降下法（ Steepest Descent Method ）

　　先に述べた基本アルゴリズムにおいて，

歩み幅： α(k) = 一定，探索方向： p(k) = -∇f(x(k)) = -[∂f/∂x₁，･･･，∂f/∂x_n]^T

とした方法です．アルゴリズムは簡単であり，ある条件の下で，大域的収束性（任意の初期点から出発しても，特定の性質を持つ点に収束）が保証されますが，収束が遅いという欠点があります．この方法の変形として，歩み幅 α(k) を 1 次元最適化によって求める方法を，最適勾配法( Optimul Gradient Method )と呼びます．

001	/************************************************/
002	/* 最急降下法による最小値の計算                 */
003	/*      f(x, y) = 100(y - x**2)**2 + (1 - x)**2 */
004	/*           最小値：(1,1) で 0.0               */
005	/*      coded by Y.Suganuma                     */
006	/************************************************/
007	#include <stdio.h>
008	#include <math.h>
009	
010	/****************************************/
011	/* xにおける微係数dxを計算する          */
012	/* （最小値のため，符号を逆にしてある） */
013	/****************************************/
014	void dsnx(double *x, double *dx)
015	{
016		dx[0] = 400.0 * x[0] * (x[1] - x[0] * x[0]) + 2.0 * (1.0 - x[0]);
017		dx[1] = -200.0 * (x[1] - x[0] * x[0]);
018	}
019	
020	/*********************************************/
021	/* 与えられた点xから，dx方向へk*dxだけ進んだ */
022	/* 点における関数値を計算する                */
023	/*********************************************/
024	double snx(double k, double *x, double *dx)
025	{
026		double w[2];
027	
028		w[0] = x[0] + k * dx[0];
029		w[1] = x[1] + k * dx[1];
030		double x1 = w[1] - w[0] * w[0];
031		double y1 = 1.0 - w[0];
032		double y  = 100.0 * x1 * x1 + y1 * y1;
033	
034		return y;
035	}
036	
037	/****************************************************************/
038	/* 黄金分割法(与えられた方向での最小値)                         */
039	/*      a,b : 初期区間 a < b                                    */
040	/*      eps : 許容誤差                                          */
041	/*      val : 関数値                                            */
042	/*      ind : 計算状況                                          */
043	/*              >= 0 : 正常終了(収束回数)                       */
044	/*              = -1 : 収束せず                                 */
045	/*      max : 最大試行回数                                      */
046	/*      w : 位置                                                */
047	/*      dw : 傾きの成分                                         */
048	/*      fun : 関数値を計算する関数の名前                        */
049	/*      return : 結果(w＋y＊dwのy)                              */
050	/****************************************************************/
051	double gold(double a, double b, double eps, double *val, int *ind, int max,
052	            double *w, double *dw, double (*fun)(double, double *, double *))
053	{
054						// 初期設定
055		*ind       = -1;
056		double x   = 0.0;
057		double tau = (sqrt(5.0) - 1.0) / 2.0;
058		double x1  = b - tau * (b - a);
059		double x2  = a + tau * (b - a);
060		double f1  = fun(x1, w, dw);
061		double f2  = fun(x2, w, dw);
062		int count  = 0;
063						// 計算
064		while (count < max && *ind < 0) {
065			count += 1;
066			if (f2 > f1) {
067				if (fabs(b-a) < eps) {
068					*ind = 0;
069					x    = x1;
070					*val = f1;
071				}
072				else {
073					b  = x2;
074					x2 = x1;
075					x1 = a + (1.0 - tau) * (b - a);
076					f2 = f1;
077					f1 = fun(x1, w, dw);
078				}
079			}
080			else {
081				if (fabs(b-a) < eps) {
082					*ind = 0;
083					x    = x2;
084					*val = f2;
085					f1   = f2;
086				}
087				else {
088					a  = x1;
089					x1 = x2;
090					x2 = b - (1.0 - tau) * (b - a);
091					f1 = f2;
092					f2 = fun(x2, w, dw);
093				}
094			}
095		}
096						// 収束した場合の処理
097		if (*ind == 0) {
098			*ind      = count;
099			double fa = fun(a, w, dw);
100			double fb = fun(b, w, dw);
101			if (fb < fa) {
102				a  = b;
103				fa = fb;
104			}
105			if (fa < f1) {
106				x    = a;
107				*val = fa;
108			}
109		}
110	
111		return x;
112	}
113	
114	/********************************************************/
115	/* 最急降下法(steepest descent method)                  */
116	/*      opt_1 : =0 : １次元最適化を行わない             */
117	/*              =1 : １次元最適化を行う                 */
118	/*      max : 最大繰り返し回数                          */
119	/*      n : 次元                                        */
120	/*      eps : 収束判定条件                              */
121	/*      step : きざみ幅                                 */
122	/*      y : 最小値                                      */
123	/*      x : x(初期値と答え)                             */
124	/*      dx : 関数の微分値                               */
125	/*      fun : 関数値を計算する関数名                    */
126	/*      dfun : 関数の微分を計算する関数名(符号を変える) */
127	/*      return : >=0 : 正常終了(収束回数)               */
128	/*               =-1 : 収束せず                         */
129	/*               =-2 : １次元最適化の区間が求まらない   */
130	/*               =-3 : 黄金分割法が失敗                 */
131	/********************************************************/
132	int steep(int opt_1, int max, int n, double eps, double step, double *y,
133	          double *x, double *dx, double (*fun)(double, double *, double *), 
134	          void (*dfun)(double *, double *))
135	{
136		double f1, f2, k, sp, y1, y2;
137		int count = 0, sw = 0;
138	
139		y1 = fun(0.0, x, dx);
140		while (count < max && sw == 0) {
141						// 傾きの計算
142			dfun(x, dx);
143			count++;
144						// １次元最適化を行わない
145			if (opt_1 == 0) {
146							// 新しい点
147				for (int i1 = 0; i1 < n; i1++)
148					x[i1] += step * dx[i1];
149							// 新しい関数値
150				y2 = fun(0.0, x, dx);
151								// 関数値の変化が大きい
152				if (fabs(y2-y1) > eps)
153					y1 = y2;
154								// 収束（関数値の変化＜ｅｐｓ）
155				else {
156					sw = count;
157					*y = y2;
158				}
159			}
160						// １次元最適化を行う
161			else {
162							// 区間を決める
163				int sw1 = 0;
164				f1 = y1;
165				sp = step;
166				f2 = fun(sp, x, dx);
167				if (f2 > f1)
168					sp = -step;
169				for (int i1 = 0; i1 < max && sw1 == 0; i1++) {
170					f2 = fun(sp, x, dx);
171					if (f2 > f1)
172						sw1 = 1;
173					else {
174						sp *= 2.0;
175						f1  = f2;
176					}
177				}
178							// 区間が求まらない
179				if (sw1 == 0)
180					sw = -2;
181							// 区間が求まった
182				else {
183								// 黄金分割法
184					k = gold(0.0, sp, eps, &y2, &sw1, max, x, dx, fun);
185								// 黄金分割法が失敗
186					if (sw1 < 0)
187						sw = -3;
188								// 黄金分割法が成功
189					else {
190									// 新しい点
191						for (int i1 = 0; i1 < n; i1++)
192							x[i1] += k * dx[i1];
193										// 関数値の変化が大きい
194						if (fabs(y1-y2) > eps)
195							y1 = y2;
196										// 収束（関数値の変化＜ｅｐｓ）
197						else {
198							sw = count;
199							*y = y2;
200						}
201					}
202				}
203			}
204		}
205	
206		if (sw == 0)
207			sw = -1;
208	
209		return sw;
210	}
211	
212	/****************/
213	/* main program */
214	/****************/
215	int main()
216	{
217				// 初期設定
218		int n       = 2;   // 変数の数
219		int max     = 10000;   // 最大試行回数
220		double eps  = 1.0e-10;   // 許容誤差
221		double step = 0.002;   // 刻み幅
222		double *x   = new double [n];
223		double *dx  = new double [n];   // 微係数
224		for (int i1 = 0; i1 < n; i1++)   // 初期位置
225			x[i1] = 0.0;
226				// 実行と出力
227		double y;
228		int sw;
229		for (int opt_1 = 0; opt_1 < 2; opt_1++) {
230					// 実行
231			sw = steep(opt_1, max, n, eps, step, &y, x, dx, snx, dsnx);
232					// 結果の出力
233			if (opt_1 == 0)
234				printf("1 次元最適化を行わなかった場合（最急降下法）\n");
235			else
236				printf("1 次元最適化を行った場合（最適勾配法）\n");
237			if (sw < 0) {
238				printf("   収束しませんでした！");
239				switch (sw) {
240					case -1:
241						printf("（収束回数）\n");
242						break;
243					case -2:
244						printf("（１次元最適化の区間）\n");
245						break;
246					case -3:
247						printf("（黄金分割法）\n");
248						break;
249				}
250			}
251			else {
252				printf("   結果＝");
253				for (int i1 = 0; i1 < n; i1++)
254					printf("%f ", x[i1]);
255				printf(" 最小値＝%f  回数＝%d\n", y, sw);
256			}
257		}
258	
259		return 0;
260	}

（出力）

1 次元最適化を行わなかった場合（最急降下法）
   結果＝0.999750 0.999500  最小値＝0.000000  回数＝8880
1 次元最適化を行った場合（最適勾配法）
   結果＝1.000000 1.000000  最小値＝0.000000  回数＝2

シンプレックス法

　　n 次元空間のシンプレックス（単体）を利用し，最適解を求める方法です．シンプレックスとは，n 次元空間において 0 でない体積を囲む少なくとも ( n + 1 ) 個の点から構成されます．例えば，4 次元空間（ 3 変数関数）においては 4 面体，3 次元空間（ 2 変数関数）においては 3 角形に相当します．シンプレックス法では，初期シンプレックスを与えた後，そのシンプレックスを移動，縮小，又は，拡大しながら最適解に収束させようとします．

　　概念的なことを理解してもらうために，2 変数関数

z = f(x, y)

の最小値を求める場合について，アルゴリズムの細部は省略して説明します．まず，3 点 A，B，C （初期シンプレックス）を決めます．次に，3 点の内，z の値が最大になる点を選び，その点を他の 2 点から決まる直線に対して対称な位置に移動します．右図においては，点 A が点 D に移動し，新しい三角形 BCD ができます．同じような処理を繰り返すことによって，

B → E，　C → F，　D → G

のように移動していきます．徐々に，Z が最小となる方向に三角形が移動していくことが分かるかと思います．しかし，このまま続けると，点 G が点 D に移動し，同じ状態を繰り返すことになってしまいますので，そのようなときは，三角形の辺の長さを半分にし，同様の処理を繰り返します．三角形の辺の長さが精度的に適切な値になったときアルゴリズムは終了します．具体的なアルゴリズムは以下に示す通りです．

n 変数の場合であれば，(n + 1) 個の初期点を選択します．一般には，1 つの初期点 x₀ を選択した後，その点の座標に，各軸方向の単位ベクトル e_i に初期値決定係数 k を乗じたものを加えた点を残りの初期点とします．
```
x_i = x₀ + k * e_i     i = 1, ･･･, n
				
```

(n + 1) 個の点から，以下に示す点を選択します．．

H : 関数値が最大になる点　　f_H = f( x_H )
G : 2 番目に関数値が最大になる点　　f_G = f( x_G )
L : 関数値が最小になる点　　f_L = f( x_L )

点 H を除いた残りの点の重心 x_C を計算し，この重心に関する点 H の鏡映，
```
x_R = 2x_C - x_H
				
```
をもとめ，その関数値 f_R = f(x_R) を計算します．
1. もし，この点がより良い点でなければ，つまり，
```
f_R ≧ f_H
					
```
  であれば，ステップが大きすぎることになるので縮小を行います．
```
x_N = ( 1 - λ₁)x_H + λ₁x_R, f_N = f( x_N )    ただし，0 < λ₁ < 1,  λ₁ ≠ 0.5
					
```
2. もし，この点がより良い点であり，かつ，
```
f_R < ( f_L + (λ₂ - 1 )f_H) / λ₂   ただし，1 < λ₂
					
```
  であれば，この方向にもっと良い点がある可能性があるため，拡張，
```
x_E = λ₂x_R - (λ₂ - 1 )x_H, f_E = f( x_E )
					
```
  を行い，
```
f_E ≦ f_R
					
```
  であれば，x_N = x_E とします．
3. 以上述べた以外の場合は，すべて，x_N = x_R とします．
もし，
```
f_N ≧ f_G
				
```
である場合は，次の探索において，前と同じことを繰り返すことになるので，シンプレックス全体を縮小します．
```
x_i = 0.5 * (x_i + x_L)   i = 1, ･･･, n+1
				
```
次の条件が満たされれば停止し，そうでなければ，ⅱ に戻ります．

001	/************************************************/
002	/* シンプレックス法による最小値の計算           */
003	/*      f(x, y) = 100(y - x**2)**2 + (1 - x)**2 */
004	/*           最小値：(1,1) で 0.0               */
005	/*      coded by Y.Suganuma                     */
006	/************************************************/
007	#include <stdio.h>
008	
009	/*******************/
010	/* 関数値の計算    */
011	/*      y = f(x)   */
012	/*      return : y */
013	/*******************/
014	double snx(double *x)
015	{
016		double x1 = x[1] - x[0] * x[0];
017		double y1 = 1.0 - x[0];
018		double y  = 100.0 * x1 * x1 + y1 * y1;
019		return y;
020	}
021	
022	/******************************************/
023	/* シンプレックス法                       */
024	/*      max : 最大繰り返し回数            */
025	/*      n : 次元（変数の数）              */
026	/*      k : 初期値設定係数                */
027	/*      r1 : 縮小比率                     */
028	/*      r2 : 拡大比率                     */
029	/*      y : 最小値                        */
030	/*      x : x(初期値と答え)               */
031	/*      fun : 関数値を計算する関数名      */
032	/*      return : >=0 : 正常終了(収束回数) */
033	/*               =-1 : 収束せず           */
034	/******************************************/
035	int simplex(int max, int n, double k, double eps, double r1, double r2, double *y,
036	            double *x, double (*fun)(double *))
037	{
038				// 初期値の設定
039		double *xg = new double [n];   // 最大点以外の重心
040		double *xr = new double [n];
041		double *xn = new double [n];
042						// シンプレックス
043		double *yy  = new double [n+1];
044		double **xx = new double * [n+1];
045		for (int i1 = 0; i1 < n+1; i1++)
046			xx[i1] = new double [n];
047		for (int i1 = 0; i1 < n+1; i1++) {
048			for (int i2 = 0; i2 < n; i2++)
049				xx[i1][i2] = x[i2];
050			if (i1 > 0)
051				xx[i1][i1-1] += k;
052			yy[i1] = fun(xx[i1]);   // シンプレックスの各点の関数値
053		}
054						// 最大値，最小値となる点
055		int fh = -1;   // 最大値となる点
056		int fl = -1;   // 最小値となる点
057		for (int i1 = 0; i1 < n+1; i1++) {
058			if (fh < 0 || yy[i1] > yy[fh])
059				fh = i1;
060			if (fl < 0 || yy[i1] < yy[fl])
061				fl = i1;
062		}
063						// 2 番目に関数値が最大になる点
064		int fg = -1;   // 2 番目に関数値が最大になる点
065		for (int i1 = 0; i1 < n+1; i1++) {
066			if (i1 != fh && (fg < 0 || yy[i1] > yy[fg]))
067				fg = i1;
068		}
069				// 最小値の計算
070		int count = 0;   // 試行回数
071		int sw    = -1;
072		while (count < max && sw < 0) {
073			count++;
074						// 最大点以外の重心の計算
075			for (int i1 = 0; i1 < n; i1++)
076				xg[i1] = 0.0;
077			for (int i1 = 0; i1 < n+1; i1++) {
078				if (i1 != fh) {
079					for (int i2 = 0; i2 < n; i2++)
080						xg[i2] += xx[i1][i2];
081				}
082			}
083			for (int i1 = 0; i1 < n; i1++)
084				xg[i1] /= n;
085						// 最大点の置き換え
086			for (int i1 = 0; i1 < n; i1++)
087				xr[i1] = 2.0 * xg[i1] - xx[fh][i1];
088			double yr = fun(xr);
089			if (yr >= yy[fh]) {   // 縮小
090				for (int i1 = 0; i1 < n; i1++)
091					xr[i1] = (1.0 - r1) * xx[fh][i1] + r1 * xr[i1];
092				yr = fun(xr);
093			}
094			else if (yr < (yy[fl]+(r2-1.0)*yy[fh])/r2) {   // 拡張
095				for (int i1 = 0; i1 < n; i1++)
096					xn[i1] = r2 * xr[i1] - (r2 - 1.0) * xx[fh][i1];
097				double yn = fun(xn);
098				if (yn <= yr) {
099					for (int i1 = 0; i1 < n; i1++)
100						xr[i1] = xn[i1];
101					yr = yn;
102				}
103			}
104			for (int i1 = 0; i1 < n; i1++)
105				xx[fh][i1] = xr[i1];
106			yy[fh] = yr;
107						// シンプレックス全体の縮小
108			if (yy[fh] >= yy[fg]) {
109				for (int i1 = 0; i1 < n+1; i1++) {
110					if (i1 != fl) {
111						for (int i2 = 0; i2 < n; i2++)
112							xx[i1][i2] = 0.5 * (xx[i1][i2] + xx[fl][i2]);
113						yy[i1] = fun(xx[i1]);
114					}
115				}
116			}
117						// 最大値，最小値，2 番目に関数値が最大になる点の計算
118			fh = -1;
119			fg = -1;
120			fl = -1;
121			for (int i1 = 0; i1 < n+1; i1++) {
122				if (fh < 0 || yy[i1] > yy[fh])
123					fh = i1;
124				if (fl < 0 || yy[i1] < yy[fl])
125					fl = i1;
126			}
127			for (int i1 = 0; i1 < n+1; i1++) {
128				if (i1 != fh && (fg < 0 || yy[i1] > yy[fg]))
129					fg = i1;
130			}
131						// 収束判定
132			double e = 0.0;
133			for (int i1 = 0; i1 < n+1; i1++) {
134				if (i1 != fl) {
135					yr  = yy[i1] - yy[fl];
136					e  += yr * yr;
137				}
138			}
139			if (e < eps)
140				sw = 0;
141		}
142	
143		if (sw == 0) {
144			sw = count;
145			*y = yy[fl];
146			for (int i1 = 0; i1 < n; i1++)
147				x[i1] = xx[fl][i1];
148		}
149	
150		delete [] yy;
151		delete [] xg;
152		delete [] xr;
153		delete [] xn;
154		for (int i1 = 0; i1 < n+1; i1++)
155			delete [] xx[i1];
156		delete [] xx;
157	
158		return sw;
159	}
160	
161	/****************/
162	/* main program */
163	/****************/
164	int main()
165	{
166				// 初期設定
167		int n      = 2;   // 変数の数
168		int max    = 1000;   // 最大試行回数
169		double k   = 1.0;   // 初期値設定係数
170		double eps = 1.0e-20;   // 許容誤差
171		double r1  = 0.7;   // 縮小比率
172		double r2  = 1.5;   // 拡大比率
173		double *x  = new double [n];
174		for (int i1 = 0; i1 < n; i1++)   // 初期値
175			x[i1] = 0.0;
176				// 実行
177		double y;
178		int sw = simplex(max, n, k, eps, r1, r2, &y, x, snx);
179				// 結果の出力
180		printf("シンプレックス法\n");
181		if (sw < 0)
182			printf("   収束しませんでした！\n");
183		else {
184			printf("   結果＝");
185			for (int i1 = 0; i1 < n; i1++)
186				printf("%f ", x[i1]);
187			printf(" 最小値＝%f  回数＝%d\n", y, sw);
188		}
189	
190		delete [] x;
191	
192		return 0;
193	}

（出力）

シンプレックス法
   結果＝1.000001 1.000001  最小値＝0.000000  回数＝82

［演習６］この講義で示した最急降下法，及び，シンプレックス法を使用して，様々な非線形関数の最小値（最大値）を求めてみよ．