Keyword: 非線形計画問題, スパース
概要
本サンプルは非線形計画問題(スパース)を解くC言語によるサンプルプログラムです。 本サンプルは以下に示される非線形関数を最小化する解を求めて出力します。
※本サンプルはnAG Cライブラリに含まれる関数 nag_opt_nlp_sparse() のExampleコードです。本サンプル及び関数の詳細情報は nag_opt_nlp_sparse のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
オプションパラメータ
(本関数の詳細はnag_opt_nlp_sparse のマニュアルページを参照)| このデータをダウンロード |
nag_opt_nlp_sparse (e04ugc) Example Program Optional Parameters Begin e04ugc minor_iter_lim = 20 iter_lim = 30 End
- 1行目はタイトル行で読み飛ばされます。
- 3行目にe04ugcのパラメータの開始を示しており、読み飛ばされます。
- 4行目にminor_iter_limの値として"20"を指定しています。
- 5行目にiter_limの値として"30"を指定しています。
- 6行目にパラメータの終わりを示しており、読み飛ばされます。
入力データ
(本関数の詳細はnag_opt_nlp_sparse のマニュアルページを参照)| このデータをダウンロード |
nag_opt_nlp_sparse (e04ugc) Example Program Data
Values of n and m
4 6
Values of ncnln, nonln and njnln
3 4 2
Values of nnz and iobj
14 6
Columns and rows names
'Varble 1' 'Varble 2' 'Varble 3' 'Varble 4' 'NlnCon 1'
'NlnCon 2' 'NlnCon 3' 'LinCon 1' 'LinCon 2' 'Free Row'
Matrix nonzeros: value, row index, column index
1.0E+25 1 1
1.0E+25 2 1
1.0E+25 3 1
1.0 5 1
-1.0 4 1
1.0E+25 1 2
1.0E+25 2 2
1.0E+25 3 2
1.0 5 2
-1.0 4 2
3.0 6 3
-1.0 1 3
-1.0 2 4
2.0 6 4
Lower bounds
-0.55 -0.55 0.0 0.0 -894.8 -894.8 -1294.8 -0.55
-0.55 -1.0E+25
Upper bounds
0.55 0.55 1200.0 1200.0 -894.8 -894.8 -1294.8 1.0E+25
1.0E+25 1.0E+25
Initial estimate of X
0.0 0.0 0.0 0.0
- 1行目はタイトル行で読み飛ばされます。
- 4行目に変数の数(n)、一般制約の数(m)を指定しています。
- 7行目に非線形制約の数(ncnln)、非線形目的関数の数(nonln)、非線形ヤコビ変数の数(njnln)を指定しています。
- 10行目にヤコビ行列Aの非ゼロ要素の数(nnz)、目的関数の線形部分の非ゼロ要素を含む行(iobj)を指定しています。
- 13〜14行目に列と行の名前(crnames)を指定しています。
- 17〜30行目にヤコビ行列の非ゼロ要素の値(a)、行のインデックス(ha)、列のインデックス(icol)を指定しています。
- 33〜34行目に変数の下限と制約の下限(bl)を指定しています。
- 37〜38行目に変数の上限と制約の上限(bu)を指定しています。
- 41行目にxの初期推定値(xs)を指定しています。
出力結果
(本関数の詳細はnag_opt_nlp_sparse のマニュアルページを参照)| この出力例をダウンロード |
nag_opt_nlp_sparse (e04ugc) Example Program Results Final objective value = 5126.498 Optimal X = 0.119 -0.396 679.945 1026.067 A run of the same example with a warm start: -------------------------------------------- Final objective value = 5126.498 Optimal X = 0.119 -0.396 679.945 1026.067
- 6行目にオプション・パラメータ用のファイルe04ugce.optが読み込まれていることが示されています。
- 8行目にminor_iter_limに"20"が設定されていることが示されています。
- 9行目にiter_limに"30"が設定されていることが示されています。
- 11〜58行目に以下に示すプログラムのオプション引数が出力されています。
fcheck 直近の主となる反復の後にこの回数の反復ごとに現在の解が一般線形制約を満たしているかを確認するための数値テストが行われます。 expand_freq 実行可能許容値がこの値に増加します。 factor_freq 基底行列の因数分解で生じる基底変換の最大数。 scale_tol スケーリングパスの数を制御します。 minor_feas_tol 全ての変数の下限・上限がこの許容値の範囲内かどうか確認されます。 scale_opt 変数や制約のスケーリング。 minor_opt_tol QP子問題の最適化の判定に使用されます。 part_price プライシングに必要な処理を軽減します。 crash_tol Crash処理で三角基底の検索時に行列Aの非ゼロ要素を無視できるようにします。 pivot_tol Pivot要素がこの値より小さい場合無視されます。 minor_print_level 小さな反復の出力結果を制御します。"Nag_NoPrint"は何も出力されないことを意味します。 crash 行列Aのどの行や列が基底に値するかを決定し、Crash処理が何回呼び出されるかを示します。 "Nag_NoCrash"は初期基底がB=Iとなるスラック変数のみを含むことを意味します。 elastic_wt 初期の重み。 obj_deriv 目的関数の勾配の全要素が関数objfunで与えられているかどうかを示します。"Nag_TRUE"は与えられていることを意味します。 con_deriv 制約のヤコビ行列の全要素が関数confunで与えられているかどうかを示します。"Nag_TRUE"は与えられていることを意味します。 verify_glad 導関数チェックのレベル。"Nag_SimpleCheck"は簡易テスト(3回のobjfunの実行と2回のconfunの実行)のみ実行されることを意味します。 print_deriv 導関数チェックの結果の出力を制御。"Nag_D_Print"は導関数チェックの詳細が出力されることを意味します。 Start obj check at col.. 関数objfunで計算された勾配要素の導関数チェックの開始のカラム。 Stop obj check at col.. 関数objfunで計算された勾配要素の導関数チェックの終了のカラム。 Start con check at col.. 関数confunで計算されたヤコビ行列の導関数チェックの開始のカラム。 Stop con check at col.. 関数confunで計算されたヤコビ行列の導関数チェックの終了のカラム。 direction 最適化の方向。"Nag_Minimize"は最適化の方向が最小化であることを意味します。 Nonlinear Objective vars 非線形目的変数。 major_opt_tol 双対変数の最終的な精度。 f_prec 非線形関数が計算できる相対的精度の基準となる相対関数。 inf_step f(x + αp)のαがこの値に達したら反復が終了します。 max_sb superbasic変数に割り当てられた記憶域の上限。 f_diff_int 導関数の推定に使用される区間。 unbounded_obj 目的関数fの絶対値がこの値に達したら反復が終了します。 c_diff_int 勾配のより正確な推定値を得るために最適解の近くで使用されます。 major_step_lim 直線探索(linesearch)の際のxの変量の上限。 deriv_linesearch 直線探索(linesearch)が三次補間を使用するかどうかを指定するフラグ。 "Nag_TRUE"は三次補間が使用されることを意味します。 print_level 結果出力のレベル。"Nag_Soln_Iter"は最終解と各反復の結果が出力されることを意味します。 major_iter_lim 大きな反復の最大数。 linesearch_tol 各反復の検索の際の精度を制御します。 minor_iter_lim 非線形制約の一連の制約化で可能な反復の最大数。 inf_bound 問題制約の定義での無限境界。 iter_lim 小さな反復の最大数。 hess_storage ラグランジュ関数のヘッセ行列の準ニュートン近似の格納及び更新の手法を指定します。 "Nag_HessianFull"は近似ヘッセ行列が密行列として処理されることを意味します。 hess_update ヘッセ行列更新ベクトルのペアの最大数。 hess_freq この回数のBFGS(Broyden-Fletcher-Goldfarb-Shanno)更新から形成される近似ヘッセ行列をリセットします。 Nonlinear constraints 非線形制約。 major_feas_tol 非線形制約の精度。 Nonlinear Jacobian vars 非線形ヤコビ変数。 violation_limit 直線探索(linesearch)後の最大制約違反の大きさの絶対限度。 Variables 変数。 Linear constraints 線形制約。 Nonlinear variables 非線形変数。 Linear variables 線形変数。 lu_factor_tol 再因数分解時の基底因数分解の安定性及びスパース性を制御します。 lu_sing_tol 悪条件の基底行列を防ぐために使用される特異点許容値。 lu_update_tol 更新時の基底因数分解の安定性及びスパース性を制御します。 lu_den_tol 基底行列のLU分解時に使用される密度許容値。 eps (machine precision). マシン精度。 start プログラムがどのように開始するかを示しています。"Nag_Cold"は初期のCrash処理が使用されることを意味します。。 feas_exit 非線形制約について実行可能点を見つけるために実行される追加の反復が必要かどうかを示しています。 "Nag_FALSE"はこれ以上反復が実行されないことを意味します。 Names 変数の名前が提供されていることを示しています。 print_80ch 大きい反復及び小さい反復で出力される行の最大の長さ。 outfile 結果が出力されるファイル名。 - 60〜63行目に以下に示すメモリの割り当てが出力されています。
nz_coef 基底因子用にどのくらいのメモリが初めに割り当てられるかを指定します。 state オプション引数stateに必要なメモリがnAGライブラリの関数によって自動的に割り当てられたことを示しています。 Integers 整数の場合の作業用配列の初期サイズ。 lambda オプション引数lambdaに必要なメモリがnAGライブラリの関数によって自動的に割り当てられたことを示しています。 Reals 実数の場合の作業用配列の初期サイズ。 - 65〜67行目にはスケールオプションが1から0に減少したことと、実行不可能点の合計が出力されています。
- 69行目には関数confunが6つの制約勾配から6つを設定したことが出力されています。
- 70行目には関数objfunが4つの目的関数の勾配から4つを設定したことが出力されています。
- 78行目にはヤコビ行列に問題がないことが出力されています。
- 80行目には制約の最大の相違が出力されています。
- 89行目には目的関数の勾配に問題がないことが出力されています。
- 90行目には推定される勾配が出力されています。
- 91行目には差分近似が出力されています。
- 94行目には弾性モード開始時の重みが出力されています。
- 97〜107行目には以下が出力されています。
Maj 大きな反復の回数。 Mnr 実行可能解生成フェーズや最適化フェーズで必要とされる小さな反復の回数。 Step 探索方向に進んだステップ幅。 Merit Function ラグランジュメリット関数の値。 Feasibl スケーリングされた非線形制約残差ベクトルの最大要素の値。 Optimal 最大の補間ギャップベクトルの最大要素の値。 Cond Hz ラグランジュの縮小ヘッセ行列の条件数の推定値。 PD 実行可能性と最適性についての収束判定の状況。 "T"は満たしていることを意味し、"F"は満たしていないことを意味します。さらに右側に表示されている"R"は近似ヘッセ行列がリセットされたことを意味し、"i"はQP子問題が実行不可能であることを意味し、"s"は自己スケーリング BFGS(Broyden-Fletcher-Goldfarb-Shanno)更新が実行されたことを意味します。 - 109〜110行目にはNP問題を解くのに9回の大きな反復と19回の小さな反復が行われたことが示されています。
- 112〜117行目には以下が出力されています。
Variable 変数の名前。 State 変数の状態(FRは下限と上限が等しいことを意味します)。 Value 最後の反復の変数の値。 Lower Bound 下限。 Upper Bound 上限。 Lagr Mult 関連する境界制約のラグランジュ乗数。 Residual 残差(変数の値と上限/下限の近いほうとの差)。 - 119〜126行目には以下が出力されています。
Constrnt 制約の名前。 State 制約の状態(ULは上限であることを意味し、LLは下限であることを意味します)。 Value 最後の反復の制約の値。 Lower Bound 下限。 Upper Bound 上限。 Lagr Mult 関連する境界制約のラグランジュ乗数。 Residual 残差(制約の値と上限/下限の近いほうとの差)。 - 128行目に最適解が見つかったことが示されています。
- 130行目に最終的な目的関数値が出力されています。
- 134行目からはCold startと同じ例を用いたWarm startの場合の実行結果が出力されています。
- 137〜184行目にプログラムのオプション引数が出力されています。
- 186〜189行目にメモリの割り当てが出力されています。
- 197〜205行目には以下が出力されています。
Maj 大きな反復の回数。 Mnr 実行可能解生成フェーズや最適化フェーズで必要とされる小さな反復の回数。 Step 探索方向に進んだステップ幅。 Merit Function ラグランジュメリット関数の値。 Feasibl スケーリングされた非線形制約残差ベクトルの最大要素の値。 Optimal 最大の補間ギャップベクトルの最大要素の値。 Cond Hz ラグランジュの縮小ヘッセ行列の条件数の推定値。 PD 実行可能性と最適性についての収束判定の状況。 "T"は満たしていることを意味し、"F"は満たしていないことを意味します。さらに右側に表示されている"R"は近似ヘッセ行列がリセットされたことを意味し、"S"は自己スケーリング BFGS(Broyden-Fletcher-Goldfarb-Shanno)更新の修正が必要だったことを意味します。 - 207〜208行目にはNP問題を解くのに7回の大きな反復と8回の小さな反復が行われたことが示されています。
- 210行目に最適解が見つかったことが示されています。
- 212行目に最終的な目的関数値が出力されています。
ソースコード
(本関数の詳細はnag_opt_nlp_sparse のマニュアルページを参照)
※本サンプルソースコードはnAG数値計算ライブラリ(Windows, Linux, MAC等に対応)の関数を呼び出します。
サンプルのコンパイル及び実行方法
| このソースコードをダウンロード |
/* nag_opt_nlp_sparse (e04ugc) Example Program.
*
* CLL6I261D/CLL6I261DL Version.
*
* Copyright 2017 Numerical Algorithms Group.
*
* nAG C Library
*
* Mark 26.1, 2017.
*
*/
#include <nag.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <nag_stdlib.h>
#include <nage04.h>
#ifdef __cplusplus
extern "C"
{
#endif
static void nAG_CALL confun(Integer ncnln, Integer njnln,
Integer nnzjac, const double x[], double conf[],
double conjac[], Nag_Comm *comm);
static void nAG_CALL objfun(Integer nonln,
const double x[], double *objf,
double objgrad[], Nag_Comm *comm);
#ifdef __cplusplus
}
#endif
#define NAMES(I, J) names[(I)*9+J]
int main(void)
{
const char *optionsfile = "e04ugce.opt";
Integer exit_status = 0, *ha = 0, i, icol, iobj, j, jcol, *ka = 0, m, n,
ncnln;
Integer ninf, njnln, nnz, nonln;
Nag_E04_Opt options;
char **crnames = 0, *names = 0;
double *a = 0, *bl = 0, *bu = 0, obj, sinf, *xs = 0;
Integer verbose_output;
Nag_Comm comm;
NagError fail;
INIT_FAIL(fail);
printf("nag_opt_nlp_sparse (e04ugc) Example Program Results\n");
fflush(stdout);
/* Skip heading in data file */
scanf(" %*[^\n]");
/* Read the problem dimensions */
scanf(" %*[^\n]");
scanf("%ld%ld", &n, &m);
/* Read NCNLN, NONLN and NJNLN from data file. */
scanf(" %*[^\n]");
scanf("%ld%ld%ld", &ncnln, &nonln, &njnln);
/* Read NNZ, IOBJ */
scanf(" %*[^\n]");
scanf("%ld%ld", &nnz, &iobj);
if (!(a = nAG_ALLOC(nnz, double)) ||
!(bl = nAG_ALLOC(n + m, double)) ||
!(bu = nAG_ALLOC(n + m, double)) ||
!(xs = nAG_ALLOC(n + m, double)) ||
!(ha = nAG_ALLOC(nnz, Integer)) ||
!(ka = nAG_ALLOC(n + 1, Integer)) ||
!(crnames = nAG_ALLOC(n + m, char *)) ||
!(names = nAG_ALLOC((n + m) * 9, char))
)
{
printf("Allocation failure\n");
exit_status = 1;
goto END;
}
/* Read the column and row names */
scanf(" %*[^\n]");
scanf(" %*[^']");
for (i = 0; i < n + m; ++i) {
scanf(" '%8c'", &NAMES(i, 0));
NAMES(i, 8) = '\0';
crnames[i] = &NAMES(i, 0);
}
/* read the matrix and set up ka. */
jcol = 1;
ka[jcol - 1] = 0;
scanf(" %*[^\n]");
for (i = 0; i < nnz; ++i) {
/* a[i] stores (ha[i], icol) element of matrix */
scanf("%lf%ld%ld", &a[i], &ha[i], &icol);
if (icol < jcol) {
/* Elements not ordered by increasing column index. */
printf("Element in column%5ld found after element in"
" column%5ld. Problem abandoned.\n", icol, jcol);
exit_status = 1;
goto END;
}
else if (icol == jcol + 1) {
/* Index in a of the start of the icol-th column equals i. */
ka[icol - 1] = i;
jcol = icol;
}
else if (icol > jcol + 1) {
/* Index in a of the start of the icol-th column equals i,
* but columns jcol+1,jcol+2,...,icol-1 are empty. Set the
* corresponding elements of ka to i.
*/
for (j = jcol + 1; j <= icol - 1; ++j)
ka[j - 1] = i;
ka[icol - 1] = i;
jcol = icol;
}
}
ka[n] = nnz;
if (n > icol) {
/* Columns N,N-1,...,ICOL+1 are empty. Set the
* corresponding elements of ka accordingly. */
for (j = icol; j <= n - 1; ++j)
ka[j] = nnz;
}
/* Read the bounds */
scanf(" %*[^\n]");
for (i = 0; i < n + m; ++i)
scanf("%lf", &bl[i]);
scanf(" %*[^\n]");
for (i = 0; i < n + m; ++i)
scanf("%lf", &bu[i]);
/* Read the initial estimate of x */
scanf(" %*[^\n]");
for (i = 0; i < n; ++i)
scanf("%lf", &xs[i]);
/* Initialize the options structure */
/* nag_opt_init (e04xxc).
* Initialization function for option setting
*/
nag_opt_init(&options);
/* Set this to 1 to cause e04ugc to produce intermediate
progress output */
verbose_output = 0;
/* Read some option values from standard input */
/* nag_opt_read (e04xyc).
* Read options from a text file
*/
nag_opt_read("e04ugc", optionsfile, &options,
(Nag_Boolean)(verbose_output==1), "stdout", &fail);
/* Set some other options directly */
options.major_iter_lim = 100;
options.crnames = crnames;
if (verbose_output == 0)
{
options.list = Nag_FALSE;
options.print_level = Nag_NoPrint;
options.print_deriv = Nag_D_NoPrint;
}
/* Solve the problem. */
/* nag_opt_nlp_sparse (e04ugc), see above. */
nag_opt_nlp_sparse(confun, objfun, n, m,
ncnln, nonln, njnln, iobj, nnz,
a, ha, ka, bl, bu, xs,
&ninf, &sinf, &obj, &comm, &options, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_opt_nlp_sparse (e04ugc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
if (verbose_output == 0)
{
printf("\n");
printf("Final objective value = %12.3f\n", obj);
printf("Optimal X = ");
for (i = 1; i <= n; ++i) {
printf("%9.3f%s", xs[i - 1], i % 7 == 0 || i == n ? "\n" : " ");
}
}
/* We perturb the solution and solve the
* same problem again using a warm start.
*/
printf("\nA run of the same example with a warm start:\n");
printf("--------------------------------------------\n");
options.start = Nag_Warm;
/* Modify some printing options */
options.print_deriv = Nag_D_NoPrint;
if (verbose_output == 1)
options.print_level = Nag_Iter;
else
options.print_level = Nag_NoPrint;
/* Perturb xs */
for (i = 0; i < n + m; i++)
xs[i] += 0.2;
/* Reset multiplier estimates to 0.0 */
if (ncnln > 0) {
for (i = 0; i < ncnln; i++)
options.lambda[n + i] = 0.0;
}
/* Solve the problem again. */
/* nag_opt_nlp_sparse (e04ugc), see above. */
fflush(stdout);
nag_opt_nlp_sparse(confun, objfun, n, m,
ncnln, nonln, njnln, iobj, nnz,
a, ha, ka, bl, bu, xs,
&ninf, &sinf, &obj, &comm, &options, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_opt_nlp_sparse (e04ugc).\n%s\n", fail.message);
exit_status = 1;
}
if (verbose_output == 0)
{
printf("Final objective value = %12.3f\n", obj);
printf("Optimal X = ");
for (i = 1; i <= n; ++i) {
printf("%9.3f%s", xs[i - 1], i % 7 == 0 || i == n ? "\n" : " ");
}
}
/* Free memory allocated by nag_opt_nlp_sparse (e04ugc) to pointers in options
*/
/* nag_opt_free (e04xzc).
* Memory freeing function for use with option setting
*/
nag_opt_free(&options, "all", &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_opt_free (e04xzc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
END:
nAG_FREE(a);
nAG_FREE(bl);
nAG_FREE(bu);
nAG_FREE(xs);
nAG_FREE(ha);
nAG_FREE(ka);
nAG_FREE(crnames);
nAG_FREE(names);
return exit_status;
}
/* Subroutine */
static void nAG_CALL confun(Integer ncnln, Integer njnln, Integer nnzjac,
const double x[], double conf[], double conjac[],
Nag_Comm *comm)
{
#define CONJAC(I) conjac[(I) -1]
#define CONF(I) conf[(I) -1]
#define X(I) x[(I) -1]
/* Compute the nonlinear constraint functions and their Jacobian. */
if (comm->flag == 0 || comm->flag == 2) {
CONF(1) = sin(-X(1) - 0.25) * 1e3 + sin(-X(2) - 0.25) * 1e3;
CONF(2) = sin(X(1) - 0.25) * 1e3 + sin(X(1) - X(2) - 0.25) * 1e3;
CONF(3) = sin(X(2) - X(1) - 0.25) * 1e3 + sin(X(2) - 0.25) * 1e3;
}
if (comm->flag == 1 || comm->flag == 2) {
/* Nonlinear Jacobian elements for column 1.0 */
CONJAC(1) = cos(-X(1) - 0.25) * -1e3;
CONJAC(2) = cos(X(1) - 0.25) * 1e3 + cos(X(1) - X(2) - 0.25) * 1e3;
CONJAC(3) = cos(X(2) - X(1) - 0.25) * -1e3;
/* Nonlinear Jacobian elements for column 2.0 */
CONJAC(4) = cos(-X(2) - 0.25) * -1e3;
CONJAC(5) = cos(X(1) - X(2) - 0.25) * -1e3;
CONJAC(6) = cos(X(2) - X(1) - 0.25) * 1e3 + cos(X(2) - 0.25) * 1e3;
}
}
static void nAG_CALL objfun(Integer nonln, const double x[], double *objf,
double objgrad[], Nag_Comm *comm)
{
#define OBJGRAD(I) objgrad[(I) -1]
#define X(I) x[(I) -1]
/* Compute the nonlinear part of the objective function and its grad */
if (comm->flag == 0 || comm->flag == 2)
*objf = X(3) * X(3) * X(3) * 1e-6 + X(4) * X(4) * X(4) * 2e-6 / 3.0;
if (comm->flag == 1 || comm->flag == 2) {
OBJGRAD(1) = 0.0;
OBJGRAD(2) = 0.0;
OBJGRAD(3) = X(3) * X(3) * 3e-6;
OBJGRAD(4) = X(4) * X(4) * 2e-6;
}
}
