線形計画（スパース）及び凸２次計画問題- サンプルソースコード Fortran

ホーム > 製品 > nAG数値計算ライブラリ > nAG Fortranライブラリ > サンプルソースコード集 > 線形計画（スパース）及び凸２次計画問題

Keyword: 線形計画, スパース, 凸２次計画問題

概要

本サンプルは線形計画（スパース）及び凸２次計画問題を解くFortranによるサンプルプログラムです。本サンプルは以下に示される二次関数を最小化する解を求めて出力します。

線形計画や２次計画問題のデータ　

※本サンプルはnAG Fortranライブラリに含まれるルーチン e04nkf() のExampleコードです。本サンプル及びルーチンの詳細情報は e04nkf のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで

入力データ

このデータをダウンロード

E04NKF Example Program Data
 7  8             :Values of N and M
48  8  7  'C'  15 :Values of NNZ, IOBJ, NCOLH, START and NNAME
'        '  '        '  '        '  '        '  '        ' :End of NAMES
'...X1...'  '...X2...'  '...X3...'  '...X4...'  '...X5...'
'...X6...'  '...X7...'  '..ROW1..'  '..ROW2..'  '..ROW3..'
'..ROW4..'  '..ROW5..'  '..ROW6..'  '..ROW7..'  '..COST..' :End of CRNAME
    0.02   7   1
    0.02   5   1
    0.03   3   1
    1.00   1   1
    0.70   6   1
    0.02   4   1
    0.15   2   1
 -200.00   8   1
    0.06   7   2
    0.75   6   2
    0.03   5   2
    0.04   4   2
    0.05   3   2
    0.04   2   2
    1.00   1   2
-2000.00   8   2
    0.02   2   3
    1.00   1   3
    0.01   4   3
    0.08   3   3
    0.08   7   3
    0.80   6   3
-2000.00   8   3
    1.00   1   4
    0.12   7   4
    0.02   3   4
    0.02   4   4
    0.75   6   4
    0.04   2   4
-2000.00   8   4
    0.01   5   5
    0.80   6   5
    0.02   7   5
    1.00   1   5
    0.02   2   5
    0.06   3   5
    0.02   4   5
-2000.00   8   5
    1.00   1   6
    0.01   2   6
    0.01   3   6
    0.97   6   6
    0.01   7   6
  400.00   8   6
    0.97   7   7
    0.03   2   7
    1.00   1   7
  400.00   8   7                                           :End of matrix A
 0.0       0.0       4.0E+02   1.0E+02  0.0      0.0      0.0       2.0E+03 
-1.0E+25  -1.0E+25  -1.0E+25  -1.0E+25  1.5E+03  2.5E+02 -1.0E+25 :End of BL 
 2.0E+02   2.5E+03   8.0E+02   7.0E+02  1.5E+03  1.0E+25  1.0E+25   2.0E+03   
 6.0E+01   1.0E+02   4.0E+01   3.0E+01  1.0E+25  3.0E+02  1.0E+25 :End of BU
 0    0    0    0    0    0    0    0   :End of ISTATE 
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0 :End of XS

１行目はタイトル行で読み飛ばされます。
２行目に変数の数（n=7）、一般制約の数（m=8）を指定しています。
３行目に線形制約行列の非ゼロ要素の数（nnz=48）、線形目的項c^Txのベクトルcの非ゼロ要素の行（iobj=8）、ヘッセ行列の非ゼロカラムの数（ncolh=7）、初期基底がどのように得られるか（start='C':内部のCrash処理が使用される）、カラム名と行名の数（nname=15）を指定しています。
４行目に線形計画問題のMPSX形式に関連する名前(names)を指定しています。
５～７行目に列と行の名前(crname)を指定しています。
８～５５行目に線形制約行列の非ゼロ要素の値(a)、行のインデックス(ha)、列のインデックス(icol)を指定しています。
５６～５７行目に変数の下限と制約の下限(bl)を指定しています。
５８～５９行目に変数の上限と制約の上限(bu)を指定しています。
６０行目にxの初期状態(istate)を指定しています。
６１行目にxの初期推定値(xs)を指定しています。

出力結果

（本ルーチンの詳細はe04nkf のマニュアルページを参照）

この出力例をダウンロード

 E04NKF Example Program Results

 Workspace provided is                 IZ(       1),  Z(       1).
 To start solving the problem we need  IZ(     428),  Z(     358).

 Exit E04NKF - Not enough integer workspace to start solving the problem.

 *** E04NKF

 Parameters
 ----------

 Frequencies.
 Check frequency.........        60       Expand frequency........     10000
 Factorization frequency.       100

 LP Parameters.
 Scale tolerance.........  9.00E-01       Feasibility tolerance...  1.00E-06
 Iteration limit.........        75       Scale option............         2
 Optimality tolerance....  1.00E-06       Partial price...........        10
 Crash tolerance.........  1.00E-01       Pivot tolerance.........  2.04E-11
 Crash option............         2

 QP objective.
 Objective variables.....         7       Hessian columns.........         7
 Superbasics limit.......         7

 Miscellaneous.
 Variables...............         7       Linear constraints......         8
 LU factor tolerance.....  1.00E+02       LU update tolerance.....  1.00E+01
 LU singularity tolerance  2.04E-11       Monitoring file.........        -1
 EPS (machine precision).  1.11E-16       Print level.............        10
 Infinite bound size.....  1.00E+20       Infinite step size......  1.00E+20
 COLD start..............                 MINIMIZE................

 Workspace provided is                 IZ(     428),  Z(     358).
 To start solving the problem we need  IZ(     428),  Z(     358).


   Itn     Step  Ninf  Sinf/Objective  Norm rg
 Itn      0 -- Infeasible.
     0  0.0E+00     1    1.152891E+03  0.0E+00
     1  4.3E+02     0    0.000000E+00  0.0E+00
 Itn      1 -- Feasible point found (for      1 equality constraints).
     1  0.0E+00     0    0.000000E+00  0.0E+00
     1  0.0E+00     0    1.460000E+06  0.0E+00
 Itn      1 -- Feasible QP solution.
     2  8.7E-02     0    9.409959E+05  0.0E+00
     3  5.3E-01     0   -1.056552E+06  0.0E+00
     4  1.0E+00     0   -1.462190E+06  0.0E+00
     5  1.0E+00     0   -1.698092E+06  1.8E-12
     6  4.6E-02     0   -1.764906E+06  7.0E+02
     7  1.0E+00     0   -1.811946E+06  1.4E-12
     8  1.7E-02     0   -1.847325E+06  1.7E+02
     9  1.0E+00     0   -1.847785E+06  9.1E-13


 Variable State     Value       Lower Bound  Upper Bound   Lagr Mult  Residual

 ...X1...    LL    0.00000           .        200.00       2361.         .
 ...X2...    BS    349.399           .        2500.0     -1.2975E-12   349.4
 ...X3...   SBS    648.853       400.00       800.00     -5.7329E-13   151.1
 ...X4...   SBS    172.847       100.00       700.00      6.4970E-13   72.85
 ...X5...    BS    407.521           .        1500.0      9.1881E-13   407.5
 ...X6...    BS    271.356           .           None    -1.1928E-12   271.4
 ...X7...    BS    150.023           .           None    -1.4130E-12   150.0


 Constrnt State     Value       Lower Bound  Upper Bound   Lagr Mult  Residual

 ..ROW1..    EQ    2000.00       2000.0       2000.0     -1.2901E+04     .
 ..ROW2..    BS    49.2316          None      60.000           .      -10.77
 ..ROW3..    UL    100.000          None      100.00      -2325.         .
 ..ROW4..    BS    32.0719          None      40.000           .      -7.928
 ..ROW5..    BS    14.5572          None      30.000           .      -15.44
 ..ROW6..    LL    1500.00       1500.0          None     1.4455E+04     .
 ..ROW7..    LL    250.000       250.00       300.00      1.4581E+04     .
 ..COST..    BS  -2.988690E+06      None         None     -1.000     -2.9887E+06

 Exit E04NKF - Optimal QP solution found.

 Final QP objective value =    -1847785.

 Exit from QP problem after       9 iterations.

１０～３４行目に以下に示すプログラムのオプション引数が出力されています。

Check frequency	直近の主となる反復の後にこの回数の反復ごとに現在の解が一般線形制約を満たしているかを確認するための数値テストが行われます。
Expand frequency	この数の反復数を超えると実行可能許容値がFeasibility toleranceの値に増加します。
Factorization frequency	基底行列の因数分解で生じる基底変換の最大数。
Scale tolerance	スケーリングパスの数を制御します。
Feasibility tolerance	実行可能点での各制約の許容可能な最大の絶対的違反。
Iteration limit	最大反復数。
Scale option	変数と制約のスケーリング。"2"は付加的なスケーリングが実行されることを意味しています。
Optimality tolerance	縮小勾配の大きさの判断に使用されます。
Partial price	プライシングに必要な処理を軽減します。
Crash tolerance	Crash処理で三角基底の検索時に行列Aの非ゼロ要素を無視できるようにします。
Pivot tolerance	Pivot要素がこの値より小さい場合無視されます。
Crash option	制約行列Aのどの行や列が基底に値するかを決定し、Crash処理が何回呼び出されるかを示します。
Objective variables	目的関数の数。
Hessian columns	ヘッセ行列のカラム数。
Superbasics limit	QP問題の非線形の度合い。
Variables	変数の数。
Linear constraints	線形制約の数。
LU factor tolerance	再因数分解時の基底因数分解の安定性及びスパース性を制御します。
LU update tolerance	更新時の基底因数分解の安定性及びスパース性を制御します。
LU singularity tolerance	悪条件の基底行列を防ぐために使用される特異点許容値。
Monitoring file	モニタリング情報が出力される論理装置番号。"-1"はモニタリング情報を出力しないことを意味しています。
EPS (machine precision)	マシン精度。
Print level	出力レベル。"10"は最終解と各反復のサマリーが出力されることを意味しています。
Infinite bound size	問題の制約における無限境界。
Infinite step size	無限解へのステップと見なされる変数の変化の大きさ。
COLD start	プログラムがどのように開始するかを示しています。"COLD start"は初期のCrash処理が使用されることを意味します。
MINIMIZE	最適化の方向が最小化であることを意味します。

３６～３７行目には与えられた作業領域の大きさと問題を解くのに必要な作業領域の大きさが出力されています。

４０～５５行目には以下が出力されています。

Itn	反復回数。
Step	探索方向に進んだステップの長さ。
Ninf	違反した制約の数。
Sinf/Objective	xが実行可能でない場合は制約違反の大きさの合計。実行可能な場合は目的関数の値。
Norm rg	縮小勾配のユークリッドノルム。

５８～６６行目には以下が出力されています。

Variable	変数の名前。
State	変数の状態（LLは下限での非基底、SBSはsuperbasic、BSは基底であることを意味します）。
Value	最後の反復の変数の値。
Lower Bound	下限。
Upper Bound	上限。
Lagr Mult	関連する境界のラグランジュ乗数。
Residual	残差（変数の値と上限／下限の近いほうとの差）。

６９～７８行目には以下が出力されています。

Constrnt	線形制約の名前。
State	線形制約の状態（ULは上限での非基底、LLは下限での非基底、EQは非基底で固定、BSは基底であることを意味します）。
Value	最後の反復の制約の値。
Lower Bound	下限。
Upper Bound	上限。
Lagr Mult	関連する境界のラグランジュ乗数。
Residual	残差（制約の値と上限／下限の近いほうとの差）。

８０行目に最適解が見つかったことが示されています。
８２行目に最終的な目的関数値が出力されています。
８４行目にはQP問題を解くのに９回の反復が行われたことが示されています。

ソースコード

（本ルーチンの詳細はe04nkf のマニュアルページを参照）

※本サンプルソースコードは科学技術・統計計算ライブラリである「nAG Fortranライブラリ」のルーチンを呼び出します。
サンプルのコンパイル及び実行方法

このソースコードをダウンロード

!   E04NKF Example Program Text
!   Mark 23 Release. nAG Copyright 2011.

    MODULE e04nkfe_mod

!      E04NKF 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 qphx(nstate,ncolh,x,hx)

!         Routine to compute H*x. (In this version of QPHX, the Hessian
!         matrix H is not referenced explicitly.)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          INTEGER, INTENT (IN)                :: ncolh, nstate
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: hx(ncolh)
          REAL (KIND=nag_wp), INTENT (IN)     :: x(ncolh)
!         .. Executable Statements ..
          hx(1) = 2.0E0_nag_wp*x(1)
          hx(2) = 2.0E0_nag_wp*x(2)
          hx(3) = 2.0E0_nag_wp*(x(3)+x(4))
          hx(4) = hx(3)
          hx(5) = 2.0E0_nag_wp*x(5)
          hx(6) = 2.0E0_nag_wp*(x(6)+x(7))
          hx(7) = hx(6)

          RETURN

       END SUBROUTINE qphx
    END MODULE e04nkfe_mod
    PROGRAM e04nkfe

!      E04NKF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : e04nkf, nag_wp
       USE e04nkfe_mod, ONLY : nin, nout, qphx
!      .. 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,  &
                                              ncolh, ninf, nname, nnz, ns
       CHARACTER (1)                       :: start
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE     :: a(:), bl(:), bu(:), clamda(:),   &
                                              xs(:), z(:)
       INTEGER, ALLOCATABLE                :: ha(:), istate(:), iz(:), ka(:)
       CHARACTER (8), ALLOCATABLE          :: crname(:)
       CHARACTER (8)                       :: names(5)
!      .. Executable Statements ..
       WRITE (nout,*) 'E04NKF Example Program Results'
       FLUSH (nout)

!      Skip heading in data file.
       READ (nin,*)

       READ (nin,*) n, m
       READ (nin,*) nnz, iobj, ncolh, 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),crname(nname))

       READ (nin,*) names(1:5)
       READ (nin,*) crname(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)

!      Solve the QP problem.
!      First call is a workspace query

       leniz = 1
       lenz = 1
       ALLOCATE (iz(leniz),z(lenz))

       ifail = 1
       CALL e04nkf(n,m,nnz,iobj,ncolh,qphx,a,ha,ka,bl,bu,start,names,nname, &
          crname,ns,xs,istate,miniz,minz,ninf,sinf,obj,clamda,iz,leniz,z,lenz, &
          ifail)

       IF (ifail/=0 .AND. ifail/=12 .AND. ifail/=13) THEN
          WRITE (nout,99998) 'Query call to E04NKF failed with IFAIL =', ifail
          GO TO 20
       END IF

       DEALLOCATE (iz,z)

       lenz = minz
       leniz = miniz
       ALLOCATE (iz(leniz),z(lenz))

       ifail = 0
       CALL e04nkf(n,m,nnz,iobj,ncolh,qphx,a,ha,ka,bl,bu,start,names,nname, &
          crname,ns,xs,istate,miniz,minz,ninf,sinf,obj,clamda,iz,leniz,z,lenz, &
          ifail)

20     CONTINUE

99999  FORMAT (/1X,A,I5,A,I5,A,A)
99998  FORMAT (1X,A,I5)
    END PROGRAM e04nkfe

線形計画（スパース）及び凸２次計画問題

Fortranによるサンプルソースコード : 使用ルーチン名：e04nkf

概要

入力データ

出力結果

ソースコード