森林レベルの最適化
吉本 敦
統計数理研究所

数理計画法による最適化モデル

線形計画法LP
Strategic Planning

0-1 整数計画法
空間制約・場所特定
Tactical Planning0. オーソドックスな問題：持続性（時間）複数林分の管理
制御時期，場所，強度：時空間を考慮
1415191617202223211821061171281391345

数理計画法の枠組みで定式化• 決定変数を決める：何を制御するのか？
• 必要な係数を推定する（予測値など）• 最適化の目的を決める• 制約条件を列挙する– 土地利用に関わる制約– 生産量に関わる制約：持続性– など

Production

伝統的な持続性の確保＝＞伐採量の維持Time(100 yr)

同種
同齢

同じ成長パターン

1415191617202221210611712813913451823

Strata Approach：問題の単純化

決定変数: 施業
各期にどれだけの伐採
x i,5x i, 4x i,3x i,2x i,1

-

---

X---

X---

X---

X---

X-

--

-

1st

2nd

3rd

4th

5th

Period∑x ≤ 1ji,j各分類でいつ，どれだけ？
1st

2nd

3rd

4th

5th

10%

30%

10%20%

30%

5%

25%

35%

25%

10%

30%

0%

30%

30%

10%

50%

5%

10%

10%

25%

10%

20%

30%

40%

0%

Time

Definition of variables & coefficients
Decision variable matrix
withSingle
Harvest⎛ x⎞⎟… x⎜⎜ 1,1
1,T ⎟⎟⎟⎜⎜X = ⎜    ⎟⎟,⎟⎟⎜⎜⎜⎜ x m,1  x m,T ⎟⎟⎝⎠x i,j : portion of the i-th unit to be harvested at period j

Coefficient matrix⎛ c⎞⎜⎜ 1,1 … c1,T ⎟⎟⎟⎟⎜⎜
C = ⎜    ⎟⎟⎟⎟⎜⎜⎜⎜ cm,1  cm,T ⎟⎟⎝⎠Volume flow matrix⎛⎜⎜ v1,1 … v1,T⎜V = ⎜⎜   ⎜⎜⎜⎜ v⎝ 1,m  vm,Tperiod j

⎞⎟
⎟⎟⎟⎟
⎟⎟⎟⎟
⎟⎟⎠

ci, j : coefficient of x i, j

Column vectors of all are for the same periodvi, j : volume flow from x i.j at period jObjective Function with Single Harvest
To maximize net present value of the sum of all returns fromharvests

Z = max tr(C ′X)
X⎛ c x + +c x
⎜⎜ 1,1 1,1m,1 m,1⎜⎜= max tr ⎜⎜⎜X⎜⎜
 c1,T x1,T +  + cm,T x m,T⎝m

T= max ∑ ∑ ci, j ⋅ x i, j
Xi =1 j =1

Note: trace of a square matrix is the summation of diagonal elements⎛ a⎞⎜⎜ 1,1 … a1,n ⎟⎟⎟⎟⎜⎜A = ⎜    ⎟⎟⎟⎟⎜⎜⎜⎜ an,1  an,n ⎟⎟⎝⎠⇒ntr(A) = ∑ ai,ii =1⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠

Land Accounting Constraints
with Single Harvest
Sum of portion <=0x i,1 + x i,2 x i,T −1 + x i,T ≤ 1or⎛ x⎜⎜ 1,1 … x1,T⎜⎜⎜⎜  ⎜⎜⎜⎝ x m,1  x m,Ti = 1,…,m⎞⎟⎛ ⎞⎛ 1 ⎞⎟⎟⎟⎜ 1 ⎟⎟⎜⎜ ⎟⎟⎟⎜⎜ ⎟⎟⎟⎟⎟⎜⎜  ⎟⎟ ≤ ⎜⎜⎜  ⎟⎟⎟ ⇔⎟⎟⎜⎜ ⎟⎟⎜ ⎟⎟⎟⎠⎝ 1 ⎠T ⎜⎝ 1 ⎟⎠m

X1T ≤ 1m

Harvest Flow Constraints
with Single Harvest
Even flow over timev1,j−1x1,j−1 + + vm,j−1x m,j−1 = v1,j x1,j + + vm,j x m,jor(v⎛x⎜)⎜⎜⎜ ⎜⎝ x

1, j −11, j −1

 v

m , j −1

m , j −1⎞⎟⎟⎟ = ( v⎟⎟⎠⎛ x ⎞⎟⎜)⎜⎜⎜  ⎟⎟⎟⎜⎝ x ⎟⎠
1, j1, j v

orv ′(j − 1)x (j − 1) = v ′(j)x (j)m ,jm ,j

線形計画法による定式化
Max PNVmT

Z = max tr(C ′X) = ∑ ∑ ci,j ⋅ x i,jX

最適化ソフトウェア使用
Cplex, Gurobi, SCIP(Academic Free)i=1 j =1st.
X1T ≤ 1m

土地利用制約v ′(p−1)x (p−1) = v ′(p)x (p)where,

p = 2,,T

生産量制約

v (p) = (v1,p ,v 2,p ,vm,p )′ , x (p) = (x 1,p ,x 2,p ,x m,p )′v (p)and x (p) are the p-th column vector of V and X, respectively
Note: trace of a square matrix is the summation of diagonal elements⎛ a⎞⎜⎜ 1,1 … a1,n ⎟⎟⎟⎟⎜⎜A = ⎜    ⎟⎟⎟⎟⎜⎜⎜⎜ an,1  an,n ⎟⎟⎝⎠⇒ntr(A) = ∑ ai,ii =1

複数回の伐採を考慮: Model I

Minimum Rotation: 6 periods over 10 periods

Table 1: Example of treatments

PeriodTreatment Decision

Treatment

1Coefficient

12345

678910x i,1ci,1X

0000000002x i,2ci,2X

00000X000

3x i,3ci,3X000000X00

4x i,4ci,4X

0000000X

0

5x i,5ci,5X

00000000

X

6x i,6ci,60

X

0

0

00

0000

7x i,7ci,70

X

00000

X

0

0

8x i,8ci,80

X

0

0

0000

X

0

9x i,9ci,90

X

0000000X

10xi,10ci,1000X

0000000

11x i,11ci,1100X00000X012xi,12ci,1200X000000X13xi,13ci,13000X00000014xi,14ci,14000X

00000X

15xi,15ci,15

0000

X

0000016xi,16ci,1600000X

000017xi,17ci,17000000X

00018xi,18ci,180000000X

0019xi,19ci,1900000000X020xi,20ci,20000000000XNo.

Variable

1Note: X denotes harvesting while 0 denotes no harvestingcci,j ==PNV
PNVi,jii

Growth Prediction

V3Volume (m3/ha)

300

V4V5 V6V7V2V1200

1000103050Age (year)

7090

The i-th Coefficient: ci,jPNVci,ji PNVPNVi,ji

12345678910P•V1 P•V2 P•V3 P•V4 P•V5 P•V6 P•V7 P•V8 P•V9 P•V1010PNVi, j = ∑t=1

P ⋅Vi,t ⋅ Ai
(1+ r)d⋅(t−1)PNV from the i-th variable＝PNVi,j・xi,jVolume at each period

Treatment
施業案

Periods
計画期間

11 V1
2 V1
3 V1
4 V1
5 V1
6 0
7 0
8 0
9 0
10 0
11 0
12 0
13 0
14 0
15 016 0
17 018 0
19 0
20 0
Vol Vol1

200000V2V2V2V200000000000Vol2

3000000000V3V3V300000000Vol34000000000000V4V4000000Vol4500000000000000V500000Vol56000000000000000V60000Vol670V100000000000000V7000Vol7800V2000V10000000000V800Vol8

9000V3000V200V10000000V90Vol9100000V4000V300V20600000V10Vol10

Volume flow
10Vol j = V j ⋅ ∑ Ai ⋅ xi, j

j = 1,,,10i=1

Even-Flow

Volume (106m3)

300

j-th period
200

10000246Period (years)810

Introduction of Treatment
Decision variable matrix⎛ x⎞⎜⎜ 1,1 … x1,n ⎟⎟⎟⎟⎜⎜X = ⎜    ⎟⎟⎟⎟⎜⎜⎜⎜ x m,1  x m,n ⎟⎟⎝⎠x i, j : portion of the j -th treatment is implemented for the i -th unit
Coefficient matrix⎛ c⎞⎜⎜ 1,1 … c1,n ⎟⎟⎟⎟⎜⎜C = ⎜    ⎟⎟⎟⎟⎜⎜⎜⎜ cm,1  cm,n ⎟⎟⎝⎠ci, j : coefficient of x i, j

Volume flow matrix at period p⎛ (p)(p) ⎞⎟⎟⎜⎜ v1,1 … v1,n⎟⎟⎜⎜
Vp = ⎜    ⎟⎟vi,(p)j : volume flow from x i.j⎟⎟⎜⎜(p)(p)⎟⎟⎜⎜ vv
m,n ⎟⎝ 1,m⎠Model I: Linear Programming Formulation
Max PNVmn

Z = max tr(C ′X) = ∑ ∑ ci, j ⋅ x i, j
Xi =1 j =1st.
X1n ≤ 1m

Land Accounting Constraints
mntr(Vp′X) = ∑ ∑ vi,(p)j ⋅ x i, j = v 0 , p = 1,,Ti =1 j =1

Harvest Even-Flow Constraints空間的配置：Tactical Planning
いつ、どれだけ、どこで?
位置を考慮する必要がある
同じ成長
伐採 – 伐採なし
(0-1)

同樹種
同齢
同林分
1415191617202223211821061171281391345場所特定のスケジューリング
1st period

5th period
bafeedcfighjighbadckjbeedighfhjbakedcfighjig

3rd period

cfdcba4th perioda

2nd periodkkjk

時空間的配置0-1 整数計画法による定式化

Decision Variable
Continuous => Discrete時空間的配置0-1 整数計画法による定式化
Max PNVmn

Z = max tr(C ′X) = ∑ ∑ ci,j ⋅ x i,jX

i=1 j =1st.
X1n ≤ 1m Land Accounting Constraints(1 − α)v 0 ≤ tr(Vp′X) ≤ (1 + α)v 0,p = 1,,THarvest Flow Constraints⎧⎪1x i, j = ⎪⎨⎪⎪⎩0

if the j -th treatment is implemented for the i -th standotherwise