Keyword: 非線形計画問題, スパース
概要
本サンプルは非線形計画問題（スパース）を解くFortranによるサンプルプログラムです。 本サンプルは以下に示される非線形関数を最小化する解を求めて出力します。
※本サンプルはnAG Fortranライブラリに含まれるルーチン e04ugf() のExampleコードです。本サンプル及びルーチンの詳細情報は e04ugf のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本ルーチンの詳細はe04ugf のマニュアルページを参照）| このデータをダウンロード |
E04UGF Example Program Data
4 6 :Values of N and M
3 4 2 :Values of NCNLN, NONLN and NJNLN
14 6 'C' 10 :Values of NNZ, IOBJ, START and NNAME
'Varble 1' 'Varble 2' 'Varble 3' 'Varble 4' 'NlnCon 1'
'NlnCon 2' 'NlnCon 3' 'LinCon 1' 'LinCon 2' 'Free Row' :End of NAMES
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 :End of matrix A
-0.55 -0.55 0.0 0.0 -894.8 -894.8 -1294.8 -0.55
-0.55 -1.0E+25 :End of BL
0.55 0.55 1200.0 1200.0 -894.8 -894.8 -1294.8 1.0E+25
1.0E+25 1.0E+25 :End of BU
0 0 0 0 :End of ISTATE
0.0 0.0 0.0 0.0 :End of XS
0.0 0.0 0.0 :End of CLAMDA
- １行目はタイトル行で読み飛ばされます。
- ２行目に変数の数（n=4）、一般制約の数（m=6）を指定しています。
- ３行目に非線形制約の数（ncnln=3）、非線形目的関数の数（nonln=4）、非線形ヤコビ変数の数（njnln=2）を指定しています。
- ４行目にヤコビ行列Aの非ゼロ要素の数（nnz=14）、目的関数の線形部分の非ゼロ要素を含む行（iobj=6）、初期基底がどのように得られるか（start='C':内部のCrash処理が使用される）、列名と行名の数（nname=10）を指定しています。
- ５〜６行目に列と行の名前(names)を指定しています。
- ７〜２０行目にヤコビ行列の非ゼロ要素の値(a)、行のインデックス(ha)、列のインデックス(icol)を指定しています。
- ２１〜２２行目に変数の下限と制約の下限(bl)を指定しています。
- ２３〜２４行目に変数の上限と制約の上限(bu)を指定しています。
- ２５行目に変数xの初期状態(istate)を指定しています。
- ２６行目に変数xとスラックの初期値(xs)を指定しています。
- ２７行目に非線形制約のラグランジュ乗数(clamda)を指定しています。
出力結果
（本ルーチンの詳細はe04ugf のマニュアルページを参照）| この出力例をダウンロード |
E04UGF Example Program Results
Workspace provided is IZ( 500), Z( 500).
To start solving the problem we need IZ( 628), Z( 758).
Exit E04UGF - Not enough integer workspace to start solving the problem.
*** E04UGF
Parameters
----------
Frequencies.
Check frequency......... 60 Expand frequency....... 10000
Factorization frequency. 50
QP subproblems.
Scale tolerance......... 9.00E-01 Minor feasibility tol.. 1.05E-08
Scale option............ 1 Minor optimality tol... 1.05E-08
Partial price........... 1 Crash tolerance........ 1.00E-01
Pivot tolerance......... 2.04E-11 Minor print level...... 0
Crash option............ 0 Elastic weight......... 1.00E+02
The SQP method.
Minimize................
Nonlinear objective vars 4 Major optimality tol... 1.05E-08
Function precision...... 1.72E-13 Unbounded step size.... 1.00E+20
Superbasics limit....... 4 Forward difference int. 4.15E-07
Unbounded objective..... 1.00E+15 Central difference int. 5.56E-05
Major step limit........ 2.00E+00 Derivative linesearch..
Derivative level........ 3 Major iteration limit.. 1000
Linesearch tolerance.... 9.00E-01 Verify level........... 0
Minor iteration limit... 500 Major print level...... 10
Infinite bound size..... 1.00E+20 Iteration limit........ 10000
Hessian approximation.
Hessian full memory..... Hessian updates........ 99999999
Hessian frequency....... 99999999
Nonlinear constraints.
Nonlinear constraints... 3 Major feasibility tol.. 1.05E-08
Nonlinear Jacobian vars. 2 Violation limit........ 1.00E+01
Miscellaneous.
Variables............... 4 Linear constraints..... 3
Nonlinear variables..... 4 Linear variables....... 0
LU factor tolerance..... 5.00E+00 LU singularity tol..... 2.04E-11
LU update tolerance..... 5.00E+00 LU density tolerance... 6.00E-01
eps (machine precision). 1.11E-16 Monitoring file........ -1
COLD start.............. Infeasible exit........
Workspace provided is IZ( 1256), Z( 1516).
To start solving the problem we need IZ( 628), Z( 758).
Itn 0 -- Scale option reduced from 1 to 0.
Itn 0 -- Feasible linear rows.
Itn 0 -- Norm(x-x0) minimized. Sum of infeasibilities = 0.00E+00.
confun sets 6 out of 6 constraint gradients.
objfun sets 4 out of 4 objective gradients.
Cheap test on confun...
The Jacobian seems to be OK.
The largest discrepancy was 4.41E-08 in constraint 2.
Cheap test on objfun...
The objective gradients seem to be OK.
Gradient projected in two directions 0.00000000000E+00 0.00000000000E+00
Difference approximations 1.74111992322E-19 4.48742248252E-21
Itn 0 -- All-slack basis B = I selected.
Itn 7 -- Large multipliers.
Elastic mode started with weight = 2.0E+02.
Maj Mnr Step Merit Function Feasibl Optimal Cond Hz PD
0 12 0.0E+00 3.199952E+05 1.7E+00 8.0E-01 2.1E+06 FF R i
1 2 1.0E+00 2.463016E+05 1.2E+00 3.2E+03 4.5E+00 FF s
2 1 1.0E+00 1.001802E+04 3.3E-02 9.2E+01 4.5E+00 FF
3 1 1.0E+00 5.253418E+03 6.6E-04 2.5E+01 4.8E+00 FF
4 1 1.0E+00 5.239444E+03 2.0E-06 2.8E+01 1.0E+02 FF
5 1 1.0E+00 5.126208E+03 6.0E-04 5.9E-01 1.1E+02 FF
6 1 1.0E+00 5.126498E+03 4.7E-07 2.9E-02 1.0E+02 FF
7 1 1.0E+00 5.126498E+03 5.9E-10 1.5E-03 1.1E+02 TF
8 1 1.0E+00 5.126498E+03 1.2E-12 7.6E-09 1.1E+02 TT
Exit from NP problem after 8 major iterations,
21 minor iterations.
Variable State Value Lower Bound Upper Bound Lagr Mult Residual
Varble 1 BS 0.118876 -0.55000 0.55000 -1.2529E-07 0.4311
Varble 2 BS -0.396234 -0.55000 0.55000 1.9243E-08 0.1538
Varble 3 BS 679.945 . 1200.0 1.7001E-10 520.1
Varble 4 SBS 1026.07 . 1200.0 -2.1918E-10 173.9
Constrnt State Value Lower Bound Upper Bound Lagr Mult Residual
NlnCon 1 EQ -894.800 -894.80 -894.80 -4.387 3.3644E-09
NlnCon 2 EQ -894.800 -894.80 -894.80 -4.106 6.0049E-10
NlnCon 3 EQ -1294.80 -1294.8 -1294.8 -5.463 3.3551E-09
LinCon 1 BS -0.515110 -0.55000 None . 3.4890E-02
LinCon 2 BS 0.515110 -0.55000 None . 1.065
Free Row BS 4091.97 None None -1.000 4092.
Exit E04UGF - Optimal solution found.
Final objective value = 5126.498
- ３〜４行目に与えられた作業領域の大きさと問題を解決するのに必要な作業領域の大きさが出力されています。
- １２〜５０行目に以下に示すプログラムのオプション引数が出力されています。
Check frequency 直近の主となる反復の後にこの回数の反復ごとに現在の解が一般線形制約を満たしているかを確認するための数値テストが行われます。 Expand frequency 実行可能許容値がこの値に増加します。 Factorization frequency 基底行列の因数分解で生じる基底変換の最大数。 Scale tolerance スケーリングパスの数を制御します。 Minor feasibility tol 全ての変数の下限・上限がこの許容値の範囲内かどうか確認されます。 Scale opttion 変数や制約のスケーリング。 Minor optimality tol QP子問題の最適化の判定に使用されます。 Partial price プライシングに必要な処理を軽減します。 Crash tolerance Crash処理で三角基底の検索時に行列Aの非ゼロ要素を無視できるようにします。 Pivot tolerance Pivot要素がこの値より小さい場合無視されます。 Minor print level 小さな反復の出力結果を制御します。"0"は何も出力されないことを意味します。 Crash option 行列Aのどの行や列が基底に値するかを決定し、Crash処理が何回呼び出されるかを示します。 "0"は初期基底がB=Iとなるスラック変数のみを含むことを意味します。 Elastic weight 初期の重み。 Minimize 最適化の方向が最小化であることを意味します。 Nonlinear Objective vars 非線形目的変数。 Major optimality tol 双対変数の最終的な精度。 Funciton precision 非線形関数が計算できる相対的精度の基準となる相対関数。 Unbounded step size f(x + αp)のαがこの値に達したら反復が終了します。 Superbasics limit superbasic変数に割り当てられた記憶域の上限。 Foward diffenrace int 導関数の推定に使用される区間。 Unbounded objective 目的関数fの絶対値がこの値に達したら反復が終了します。 Central differance int 勾配のより正確な推定値を得るために最適解の近くで使用されます。 Major step limit 直線探索（linesearch）の際のxの変量の上限。 Derivative linesearch 直線探索（linesearch）が三次補間を使用するかどうかが指定されています。 "Derivative linesearch"は三次補間が使用されることを意味します。 Derivative level どの非線形関数勾配がユーザ定義ルーチンで 提供されているかどうかが指定されています。。"3"は目的関数勾配と線形ヤコビアンの全要素が提供されることを意味します。 Major iteration limit 大きな反復の最大数。 Linesearch tolerance 各反復の検索の際の精度を制御します。 Verify level 勾配要素の有限差分チェック。"0"は簡易なチェックが行われることを意味します。 Minor iteration limit 非線形制約の一連の制約化で可能な反復の最大数。 Major print level 結果出力のレベル。"10"は最終解と各反復の結果が出力されることを意味します。 Infinite bound size 問題制約の定義での無限境界。 Iteration limit 小さな反復の最大数。 Hessian full memory ラグランジュ関数のヘッセ行列の準ニュートン近似の格納及び更新の手法を指定されています。 "Hessian full memory"は近似ヘッセ行列が密行列として処理されることを意味します。 Hessian updates ヘッセ行列更新ベクトルのペアの最大数。 Hessian frequency この回数のBFGS(Broyden-Fletcher-Goldfarb-Shanno)更新から形成される近似ヘッセ行列をリセットします。 Nonlinear constraints 非線形制約の数。 Major feasibility tolerance 非線形制約の精度。 Nonlinear Jacobian vars 非線形ヤコビ変数。 Violation limit 直線探索（linesearch）後の最大制約違反の大きさの絶対限度。 Variables 変数の数。 Linear constraints 線形制約の数。 Nonlinear variables 非線形変数の数。 Linear variables 線形変数の数。 LU factor tolerance 再因数分解時の基底因数分解の安定性及びスパース性を制御します。 LU singularity tolerance 悪条件の基底行列を防ぐために使用される特異点許容値。 LU update tolerance 更新時の基底因数分解の安定性及びスパース性を制御します。 LU density tolerance 基底行列のLU分解時に使用される密度許容値。 eps (machine precision). マシン精度。 Monitoring file モニタリング情報が出力される論理装置番号。"-1"はモニタリング情報が出力されないことを意味します。 COLD start プログラムがどのように開始するかを示しています。"COLD start"は初期のCrash処理が使用されることを意味します。 Infeasible exit 非線形制約を満たさない点で終了しようとしている場合、 "Infeasible exit"は追加の反復が実行されないことを意味します。 - ５２〜５３行目に与えられた作業領域の大きさと問題を解決するのに必要な作業領域の大きさが出力されています。
- ５５〜５９行目にはスケールオプションが1から0に減少したことと、実行不可能点の合計が出力されています。
- ６１行目には関数confunが６つの制約勾配から６つを設定したことが出力されています。
- ６２行目には関数objfunが４つの目的関数の勾配から４つを設定したことが出力されています。
- ６６行目にはヤコビ行列に問題がないことが出力されています。
- ６８行目には制約の最大の相違が出力されています。
- ７２行目には目的関数の勾配に問題がないことが出力されています。
- ７３行目には推定される勾配が出力されています。
- ７４行目には差分近似が出力されています。
- ７６行目にはスラック変数の基底B=Iが選択されたことが示されています。
- ７８行目には弾性モード開始時の重みが出力されています。
- ８２〜９１行目には以下が出力されています。
Maj 大きな反復の回数。 Mnr 実行可能解生成フェーズや最適化フェーズで必要とされる小さな反復の回数。 Step 探索方向に進んだステップ幅。 Merit Function ラグランジュメリット関数の値。 Feasibl スケーリングされた非線形制約残差ベクトルの最大要素の値。 Optimal 最大の補間ギャップベクトルの最大要素の値。 Cond Hz ラグランジュの縮小ヘッセ行列の条件数の推定値。 PD 実行可能性と最適性についての収束判定の状況。 "T"は満たしていることを意味し、"F"は満たしていないことを意味します。さらに右側に表示されている"R"は近似ヘッセ行列がリセットされたことを意味し、"i"はQP子問題が実行不可能であることを意味し、"s"は自己スケーリング BFGS(Broyden-Fletcher-Goldfarb-Shanno)更新が実行されたことを意味します。 - ９３〜９４行目にはNP問題を解くのに８回の大きな反復と２１回の小さな反復が行われたことが示されています。
- ９６〜１０１行目には以下が出力されています。
Variable 変数の名前。 State 変数の状態（FRは下限と上限が等しいことを意味します）。 Value 最後の反復の変数の値。 Lower Bound 下限。 Upper Bound 上限。 Lagr Mult 関連する境界制約のラグランジュ乗数。 Residual 残差（変数の値と上限／下限の近いほうとの差）。 - １０３〜１１０行目には以下が出力されています。
Constrnt 制約の名前。 State 制約の状態（ULは上限であることを意味し、LLは下限であることを意味します）。 Value 最後の反復の制約の値。 Lower Bound 下限。 Upper Bound 上限。 Lagr Mult 関連する境界制約のラグランジュ乗数。 Residual 残差（制約の値と上限／下限の近いほうとの差）。 - １１２行目に最適解が見つかったことが示されています。
- １１４行目に最終的な目的関数値が出力されています。
ソースコード
（本ルーチンの詳細はe04ugf のマニュアルページを参照）
※本サンプルソースコードは科学技術・統計計算ライブラリである「nAG Fortranライブラリ」のルーチンを呼び出します。
サンプルのコンパイル及び実行方法
| このソースコードをダウンロード |
! E04UGF Example Program Text
! Mark 23 Release. nAG Copyright 2011.
MODULE e04ugfe_mod
! E04UGF Example Program Module:
! Parameters and User-defined Routines
! .. Use Statements ..
USE nag_library, ONLY : nag_wp
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Parameters ..
INTEGER, PARAMETER :: nin = 5, nout = 6
CONTAINS
SUBROUTINE confun(mode,ncnln,njnln,nnzjac,x,f,fjac,nstate,iuser,ruser)
! Computes the nonlinear constraint functions and their Jacobian.
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Scalar Arguments ..
INTEGER, INTENT (INOUT) :: mode
INTEGER, INTENT (IN) :: ncnln, njnln, nnzjac, nstate
! .. Array Arguments ..
REAL (KIND=nag_wp), INTENT (OUT) :: f(ncnln)
REAL (KIND=nag_wp), INTENT (INOUT) :: fjac(nnzjac), ruser(*)
REAL (KIND=nag_wp), INTENT (IN) :: x(njnln)
INTEGER, INTENT (INOUT) :: iuser(*)
! .. Intrinsic Functions ..
INTRINSIC cos, sin
! .. Executable Statements ..
IF (mode==0 .OR. mode==2) THEN
f(1) = 1000.0E+0_nag_wp*sin(-x(1)-0.25E+0_nag_wp) + &
1000.0E+0_nag_wp*sin(-x(2)-0.25E+0_nag_wp)
f(2) = 1000.0E+0_nag_wp*sin(x(1)-0.25E+0_nag_wp) + &
1000.0E+0_nag_wp*sin(x(1)-x(2)-0.25E+0_nag_wp)
f(3) = 1000.0E+0_nag_wp*sin(x(2)-x(1)-0.25E+0_nag_wp) + &
1000.0E+0_nag_wp*sin(x(2)-0.25E+0_nag_wp)
END IF
IF (mode==1 .OR. mode==2) THEN
! Nonlinear Jacobian elements for column 1.
fjac(1) = -1000.0E+0_nag_wp*cos(-x(1)-0.25E+0_nag_wp)
fjac(2) = 1000.0E+0_nag_wp*cos(x(1)-0.25E+0_nag_wp) + &
1000.0E+0_nag_wp*cos(x(1)-x(2)-0.25E+0_nag_wp)
fjac(3) = -1000.0E+0_nag_wp*cos(x(2)-x(1)-0.25E+0_nag_wp)
! Nonlinear Jacobian elements for column 2.
fjac(4) = -1000.0E+0_nag_wp*cos(-x(2)-0.25E+0_nag_wp)
fjac(5) = -1000.0E+0_nag_wp*cos(x(1)-x(2)-0.25E+0_nag_wp)
fjac(6) = 1000.0E+0_nag_wp*cos(x(2)-x(1)-0.25E+0_nag_wp) + &
1000.0E+0_nag_wp*cos(x(2)-0.25E+0_nag_wp)
END IF
RETURN
END SUBROUTINE confun
SUBROUTINE objfun(mode,nonln,x,objf,objgrd,nstate,iuser,ruser)
! Computes the nonlinear part of the objective function and its
! gradient
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Scalar Arguments ..
REAL (KIND=nag_wp), INTENT (OUT) :: objf
INTEGER, INTENT (INOUT) :: mode
INTEGER, INTENT (IN) :: nonln, nstate
! .. Array Arguments ..
REAL (KIND=nag_wp), INTENT (INOUT) :: objgrd(nonln), ruser(*)
REAL (KIND=nag_wp), INTENT (IN) :: x(nonln)
INTEGER, INTENT (INOUT) :: iuser(*)
! .. Executable Statements ..
IF (mode==0 .OR. mode==2) THEN
objf = 1.0E-6_nag_wp*x(3)**3 + 2.0E-6_nag_wp*x(4)**3/ &
3.0E+0_nag_wp
END IF
IF (mode==1 .OR. mode==2) THEN
objgrd(1) = 0.0E+0_nag_wp
objgrd(2) = 0.0E+0_nag_wp
objgrd(3) = 3.0E-6_nag_wp*x(3)**2
objgrd(4) = 2.0E-6_nag_wp*x(4)**2
END IF
RETURN
END SUBROUTINE objfun
END MODULE e04ugfe_mod
PROGRAM e04ugfe
! E04UGF Example Main Program
! .. Use Statements ..
USE nag_library, ONLY : e04ugf, nag_wp
USE e04ugfe_mod, ONLY : confun, nin, nout, objfun
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Local Scalars ..
REAL (KIND=nag_wp) :: obj, sinf
INTEGER :: i, icol, ifail, iobj, jcol, &
leniz, lenz, m, miniz, minz, n, &
ncnln, ninf, njnln, nname, nnz, &
nonln, ns
CHARACTER (1) :: start
! .. Local Arrays ..
REAL (KIND=nag_wp), ALLOCATABLE :: a(:), bl(:), bu(:), clamda(:), &
xs(:), z(:)
REAL (KIND=nag_wp) :: user(1)
INTEGER, ALLOCATABLE :: ha(:), istate(:), iz(:), ka(:)
INTEGER :: iuser(1)
CHARACTER (8), ALLOCATABLE :: names(:)
! .. Intrinsic Functions ..
INTRINSIC max
! .. Executable Statements ..
WRITE (nout,*) 'E04UGF Example Program Results'
FLUSH (nout)
! Skip heading in data file.
READ (nin,*)
READ (nin,*) n, m
READ (nin,*) ncnln, nonln, njnln
READ (nin,*) nnz, iobj, start, nname
ALLOCATE (ha(nnz),ka(n+1),istate(n+m),a(nnz),bl(n+m),bu(n+m),xs(n+m), &
clamda(n+m),names(nname))
READ (nin,*) names(1:nname)
! Read the matrix A from data file. Set up KA.
jcol = 1
ka(jcol) = 1
DO i = 1, nnz
! Element ( HA( I ), ICOL ) is stored in A( I ).
READ (nin,*) a(i), ha(i), icol
IF (icol<jcol) THEN
! Elements not ordered by increasing column index.
WRITE (nout,99999) 'Element in column', icol, &
' found after element in column', jcol, '. Problem', &
' abandoned.'
GO TO 20
ELSE IF (icol==jcol+1) THEN
! Index in A of the start of the ICOL-th column equals I.
ka(icol) = i
jcol = icol
ELSE IF (icol>jcol+1) THEN
! 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.
ka((jcol+1):icol) = i
jcol = icol
END IF
END DO
ka(n+1) = nnz + 1
! Columns N,N-1,...,ICOL+1 are empty. Set the corresponding
! elements of KA accordingly.
DO i = n, icol + 1, -1
ka(i) = ka(i+1)
END DO
READ (nin,*) bl(1:(n+m))
READ (nin,*) bu(1:(n+m))
IF (start=='C') THEN
READ (nin,*) istate(1:n)
ELSE IF (start=='W') THEN
READ (nin,*) istate(1:(n+m))
END IF
READ (nin,*) xs(1:n)
IF (ncnln>0) THEN
READ (nin,*) clamda((n+1):(n+ncnln))
END IF
! Solve the problem.
! First call is a workspace query
leniz = max(500,n+m)
lenz = 500
ALLOCATE (iz(leniz),z(lenz))
ifail = 1
CALL e04ugf(confun,objfun,n,m,ncnln,nonln,njnln,iobj,nnz,a,ha,ka,bl,bu, &
start,nname,names,ns,xs,istate,clamda,miniz,minz,ninf,sinf,obj,iz, &
leniz,z,lenz,iuser,user,ifail)
IF (ifail/=0 .AND. ifail/=15 .AND. ifail/=16) THEN
WRITE (nout,99998) 'Query call to E04UGF failed with IFAIL =', ifail
GO TO 20
END IF
DEALLOCATE (iz,z)
! The length of the workspace required for the basis factors in this
! problem is longer than the minimum returned by the query
lenz = 2*minz
leniz = 2*miniz
ALLOCATE (iz(leniz),z(lenz))
ifail = -1
CALL e04ugf(confun,objfun,n,m,ncnln,nonln,njnln,iobj,nnz,a,ha,ka,bl,bu, &
start,nname,names,ns,xs,istate,clamda,miniz,minz,ninf,sinf,obj,iz, &
leniz,z,lenz,iuser,user,ifail)
20 CONTINUE
99999 FORMAT (/1X,A,I5,A,I5,A,A)
99998 FORMAT (1X,A,I5)
END PROGRAM e04ugfe