Planning cost-effective operational forest inventories

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Planning cost-effective operational forest inventories © The Author(s) 2024. Published by Oxford University Press on behalf of The International Biometric Society Published version Karppinen, Santeri; Ene, Liviu; Engberg Sundström, Lovisa; Karvanen, Juha Karppinen, S., Ene, L., Engberg Sundström, L., & Karvanen, J. (2024). Planning cost-effective operational forest inventories. Biometrics, 80(3), Article ujae104. https://doi.org/10.1093/biomtc/ujae104 2024 Biometrics , 2024, 80(3), ujae104 https://do i.org/10.1093/b iomtc/ujae104 B iometric Pra ctice P lanning c ost-effe ctiv e operation al forest inventories Sa n te ri Ka rppine n 1 , Liviu Ene 2 , L ovi sa En g be rg Sunds tröm 2 , Juha Ka rva ne n 1 , * 1 Depa rtme n t of Ma thema tics and Statistics, University of Jyv askyl a, Semin a arinkatu 15, Jyv askyl a 40014, Finland, 2 S k o gforsk, Upps al a Scie nce Pa rk, Upps al a SE-751 83, Sweden * Corresponding author: Juha Ka rva ne n, Depa rtme n t of Ma thema tics and Statistics, University of Jyv askyl a, Semin a arinkatu 15, Jyv askyl a 40014, Finland ( juha.t.ka rva ne [email protected]). ABSTRACT We address a Bayesia n two-s tage decision pro b lem in operational fores try whe re the inner stage considers scheduling the h arv esting to fulfil l de ma nd ta rgets a nd the oute r s tage conside rs selectin g the a ccura cy of pre-h arv est inv entories th a t ar e use d to estim ate the timber v olumes of the forest tracts. The higher a ccura cy of the inve n tory e nab les be tter s che duling de ci sions but al s o imp lies hi ghe r cos ts. We focus on the outer stage, which we formulate as a maximization of the pos te rior value of the inve n tory decision under a budget cons train t. The pos te rior value depe nds on the s o lution to the inne r s tage pro b le m a nd its computation is a nalytically in tra ctable, featurin g a n NP-ha rd bina ry opt imizat ion pro b lem within a hi gh-dime nsional in tegral. In pa rticula r, the bina ry optimiza tion pr o b lem is a speci al cas e of a ge ne r aliz ed quadr atic assi gnme n t proble m. We prese n t a practical method th at solv es the oute r s tage pro b le m with a n a ppr oxima tion which combines Monte Carlo s amp ling with a gre e dy, r andomiz ed method for the binary opt imizat ion pro b lem . We deriv e inv entory de cision s for a d atas e t of 100 Sw e dish forest tra cts a cross a ran ge of inve n tory bud gets a nd es tim ate the value of the inform ation to be obtaine d. KEYWORDS: Bayesian modeling; decision making; forestry; quadratic assi gnme n t pro b lem; s cheduling; v alue of information. 1 INTRODUCTION D e cisions on h arv esting are central in forestry. The scope of decision making varies from the strategic level, where the time scale is decades, to the optimal cutting of a single s te m (Kivine n, 2004 ; Ulvd al e t al., 2023 ). In ope rational short-te rm fores t p l anning, the s ta rting poin t is a s e t of forest tracts, that is, a s e t of forest areas that a forestry company has already contracted to be cut in the near future, for instance, in the next 6 months. In this setting, for est gr owth is irr eleva n t becaus e the p l anning horizon is short. Ultim ately, the de cision of in te res t is to choos e a s chedule for the h arv esting of the tracts such that the monthly demand for timber can be met as precisely as possible. Both exc e e ding and s ubc e e ding the de ma nd ta rge ts caus e los s es to the industry. The timber volume available in the tracts is uncertain, which brings the decision pro b le m in to the realm of s tatis tics. Before making the ha rves t scheduling decision, the unce rtain ty ca n be re duc e d by carrying out a forest inventory where the timber v olume is estim ate d with a chosen a ccura cy (the inve n tory decision) using field s urv eys (Nieuw e nh uis et al., 1999 ) and remotely s en s ed d ata (Vauhkonen e t al., 2014 ; Siipilehto e t al., 2016 ). A mor e accura te inve n tory leads to better scheduling decisions but also costs more because a la rge r sa mple size is ne e de d (Gr egoir e et al., 2016 ). Although the pro b lem of choosing an in- ve n tory method in industrial wood pr ocur e me n t is conceptually described already by Ståhl ( 1994 , p. 27–28), we are not aware of any previous works with practical s o lution s. In this w ork, w e dev e lop a Baye sia n a pproach for cos t-effective ope rational (short-te rm) fores t inve n tory p l a nning. We form ulate the inve n tory decision pro b lem as an opt imizat ion problem th at inv olv es the m aximization of the pos te rior value of the inve n tory decision under a budget cons train t. The c onc ept of pos te rior value, that is, the expected value of the new data to be c olle cte d, h as a crucial role in the inve n tory decision pro b lem . The diffe re nce of the pos te rior value a nd the prior value, that is, the expe cte d value before c olle cting ne w dat a, is called value of information (cf. Eidsvik et al., 2015 ). In forestry, value of information and relate d c onc epts h av e be e n a pplie d to long -term p l anning from the perspe ctiv e of a forest owner (Ståhl, 1994 ; Duvemo et al., 2014 ; Kangas et al., 2014 ; 2015 ). Computing the pos te rior value leads to s o lving a ge ne ral - iz ed quadr atic assi gnme n t proble m (G QA P) (Lee and Ma, 2004 ) which is known to be NP-hard (Cook, 1971 ). We propos e algorithm s that emp loy gre e dy heuristics and Monte Carlo s amp ling to find su ffic iently good s o lution s for operational use. Our approach for inv entory de cision p l anning allows for taking into ac c ount the prior uncert aintie s in the timber volume, as well as the unce rtain ties in the measure me n ts ca rried out in a fores t inve n tory. The prima ry output of our method is a p l an on how to c ost-effe ctiv e ly se le ct tract inv e n tory a ccura cies unde r unce rtain ty. Our practical exa mple i l lustrates this for 100 Sw e dish forest tracts. We begin in the next section by exp l aining the inve n tory deci - sion pro b le m of in te res t in more detail and pro vidin g guidance for the rest of the pa pe r. Web Appe ndix 1 con tains a table of all notation used in this pa pe r. Re c eiv e d: Nov e mbe r 20, 2023; Revise d: July 6, 2024; Ac c epte d: Septe mbe r 5, 2024 ©The Author(s) 2024. Published by Oxford University Press on behalf of The In te rn ation al Biome tric Socie ty. Thi s i s a n Ope n Ac c ess a rticle dis tributed unde r the te rms of the C re ative Common s A ttribution Licen s e ( https://creativ ec ommons .org/lic ense s/by/4.0/ ), which permits unre stricted re use, dis tribution, a nd reproduction in any me dium, provide d the original work is properly cited. Downloaded from https://academic.oup.com/biometrics/article/80/3/ujae104/7774592 by Jyvaskyla University user on 27 September 2024 2 r Biometrics , 2024, Vol. 80, No. 3 2 OVERVIEW OF THE DECISION PROBLEM We consider a forestry company th at h as the right to cut n S forest tracts that each h av e n A timber as s ortments . The c ompany has some prior (pre-inve n tory) information on the timber volumes Vof the tracts in the form of an (n S ×n A ) -dimensional prior probability dis tribution p V . Furthe rmore, the compa ny has the opportunity to improve the kno wledg e on Vby making a decision x (I) on a forest inventory. The result of x (I) is an inventory d atas e t Y | x (I) fo llowing a n (n S ×n A ) -dime n sional pro babil ity d istribut ion p Y| x ( I) , which int egrat es ove r the prior p V a nd a measure me n t probabil ity d istribut ion p Y| V, x ( I) . The decision v e ctor x (I) encodes the selection of an inve n tory method for each tra ct amon g n I pre define d inv e n tory me thods of v aryin g a ccuracy. The cost of applying inventory method i to tract s is C s,i , a known cons ta n t. Under this s e tting, w e c onsider the c ost-effe ctiv e m aximization of the pos te rior value of the inve n tory decision x (I) which leads us to model a two-stage deci sion-mak ing pro b le m whe re a n “inne r” pro b le m is nes ted inside a n “oute r” pro b lem . The inner pro b lem con siders s che duling the h arv ests of the forest tracts in order to meet a known demand D a,t for each as s ortment a at (future) times t = 1 , 2 , . . . , n T , a nd the oute r pro b lem con siders the c ost-effe ctiv e sele ction of x (I) under a budget cons train t, B . The outer and inner pro b lem s are linke d, sinc e the s o lution of the inner pro b lem depends on the probabilistic kno wledg e available on the volumes V, which the c ompany h as the opportunity to impro ve usin g Y | x (I) obtained from the outer decision problem. Ev en though w e model the full tw o-stage pro b lem, we focus on s o lving the outer pro b lem . Thi s i s becaus e the s o lution of the inner pro b lem is dependent on which Y | x (I) is realized. In the p l anning of x (I) , Ywi l l ne c es s a rily be unknown a nd ra ndom, a nd ther efor e a s en sib le, d at a -informed s o lution for the inner proble m ca n only be obtaine d onc e (the chosen) x (I) h as be e n ca rried out in practice, and a r eal-world r ealiza tion yof the ra ndom va ri - able Y | x (I) c olle cte d. In the c omputation al method that selects the inve n tory de cision, how ev e r, the inne r pro b lem is s o lv e d for multiple r ealiza tions Y | x (I) for each candidate decision x (I) . In effect, this evaluates the fitness of a given x (I) on average, but does not provide a single, concrete solution to the inner problem. We return to this in the d isc ussion. The rest of the paper is organized as follows. Section 3 introduc es the Sw e dish tract d atas e t us ed in our practical i l lustration. The models for the prior volume distribution p V and the measure me n t dis tribut ion p Y| V, x ( I) are describe d in Se ction 4 . The inve n tory decision pro b lem is formally prese n ted as a n optimization pro b lem in Section 5 . In Se ction 6 , w e prese n t a practical method for solving the inve n tory decision pro b le m, a nd s tudy its pe rforma nce in sim ul ation s. Finally, we app ly the me thod to the Sw e dish tract d atas e t in Se ction 7 and c onclude with a discussion in Section 8 . 3 DATA The d atas e t for the practical i l lustr ation is a r a ndom sa mple of a la rge r dataset a nd e ncompasses n S = 100 forest tracts in Southern Sw e den . Thes e tracts h av e be en in fact clearcut but here we consider the pre-ha rves t situation with three inve n tory methods ( n I = 3 ), three as s ortme n ts (pine , spruce , or deciduous trees, n A = 3 ) and a 6 months’ h arv est sche duling horizon ( n T = 6 ). Letting s = 1 , 2 , . . . , n S , a = 1 , 2 , . . . , n A , i = 1 , 2 , . . . , n I , and t = 1 , 2 , . . . , n T index tracts, as s ortme n ts, inve n tory methods a nd time poin ts, respe ctiv ely, the data consist of 1. the prior means μ0 a,s and v ari ances (σ0 a,s ) 2 of timber volumes for each as s ortme n t a nd tract, 2. the measure me n t va ria nces σ2 a,s,i for the inve n tory methods i = 1 , 2 , 3 for each as s ortme n t a nd tract, 3. the costs C s,i of applying inventory method i in tract s , 4. the inve n tory bud get B that the inve n tory cos ts ca nnot exc e e d, and 5. the de ma nd ta rg ets D a,t for e ach as s ortme n t a nd mon th. The prior means and v ari ances are o btained fr om a pr eviously dev elope d imputation model, which was fitted using data on clea r-cut tracts whe re both the true timbe r volumes a nd the cha racte ris tics of the tracts are known. A more detailed description of the data and the process is given in Web Appendix 2 . The measure me n t va ria nces σ2 a,s,i of the three inve n tory methods ( i = 1 , 2 , 3 ) are als o o btained using a pre viously de v elope d model. The inve n tory methods a re cha racte rized b y the n umbe r of fixed -a r ea field plots tha t ar e use d in the estim ation of timber volume (5, 10, and 20 plots per forest tract, respe ctiv ely). The model is based on data on clear-cut tracts where the location and the volume of each tre e h av e be en re c orde d by the h arv ester. Thes e d a ta wer e us ed to simul ate inve n tories in-silico where the timber volume of the trees inside a plot is estim ate d without cutting the tre es . More details are given in Web Appendix 3 . The costs C s,i were assumed to be cons ta n t ove r the s e t of forest tracts so that C s,i = C i . We assume costs C 1 = 100 , C 2 = 150 , and C 3 = 250 that consist of a fixed cost of 50 monetary units a nd the cos t of measuring 5, 10, and 20 plots, respe ctiv ely, at a cost of 10 monetary units per plot. The inventory budget B wi l l be varied from 10 000 to 25 000 in incre me n ts of 1000. The de ma nd ta rgets a re s e t b y the indus try in the real world but in the practical i l lustration they are obtained by setting D a,t =  n S s =1 μ0 a,s /n T for each as s ortme n t a and time t. This yields for all 6 months c onsidere d the demand targets 1520.11, 7495.42, and 481.26 for pine, spruce, and deciduous trees, respe ctiv ely. 4 S TATIS TICAL M ODEL S F OR ASS ORTMENT VOLUME S A ND INVENTORY M EASUREM ENTS Next, we detail the prior distribution of the timber volumes, p V , a nd the measure me n t model p Y| V, x (I) in the description of Section 2 and d isc uss their est imat ion. In ge ne ral, we use V a,s to refer to the volume of timber as s ortme n t (eg, pine, spruce, or deciduous trees) a ∈ { 1 , 2 , . . . , n A } , in tract s ∈ { 1 , 2 , . . . , n S } , a nd de note the corresponding (random) measure me n t b y Y a,s . Furthe rmore, we wi l l use Vand Y, respe ctiv ely, t o denot e the v olumes and meas urements of all assortme n ts in all tracts (stacked to a random v e ctor in some order that p l ays no particul ar ro le). Downloaded from https://academic.oup.com/biometrics/article/80/3/ujae104/7774592 by Jyvaskyla University user on 27 September 2024 Biometrics , 2024, Vol. 80, No. 3 r 3 We model the prior (pre-inve n tory) knowled ge on the (univ ari ate) vo lume V a,s by the lo g -norm al distribution V a,s ∼log - N(ν(p) a,s , λ(p) a,s ) , (1) whe re ν(p) a,s a nd λ(p) a,s a re the pa ra mete rs of the log -norm al dis tri - bution. The distribution of the measurement Y a,s , given true volume v a,s and inventory method i , i s likewi se modeled by a log -norm al distribution: Y a,s | (V a,s = v a,s , x (I) s = i ) ∼log - N(ν(m ) a,s,i , λ(m ) a,s,i ) , (2) with pa ra mete rs ν(m ) a,s,i a nd λ(m ) a,s,i , a nd with x (I) s ∈ { 1 , 2 , . . . , n I } s ta nding for the inve n tory decision a pplied to tract s . We estimate the pa ra mete rs ν(p) a,s , λ(p) a,s , ν(m ) a,s,i , a nd λ(m ) a,s,i in ( 1 ) a nd ( 2 ) using the method of mome n ts, b y matching the mean of ( 1 ) and the v ari ances of ( 1 ) and ( 2 ) to the data values μ0 a,s , (σ0 a,s ) 2 , and σ2 a,s,i th at w ere d isc ussed in Section 3 . The details are given in Web A ppendix 4 . A pr ior i, w e ass ume th at each v olume V a,s is indepe nde n t of the volumes of other assortme n ts in the same or other tra cts. S imil arly, the o bs erv ation s Y a,s are c ondition ally independent given V a,s = v a,s and the inventory decision. We d isc uss extensions of this model in Section 8 . Finally, Web Appendix 5 also prese n ts a c omputation ally simple r model whe re the prior ( 1 ) a nd measure me n t model ( 2 ) ar e r ep l ac e d with norm al distribution s in stead. While this model ca n be unrealis tic in the se nse that it places probability mass on ne gativ e o bs erv ation s, it may sti l l be useful with moderate prior a nd measure me n t va ria nc es . 5 THE FOREST INVENTORY PLANNING PROBLEM The p l anning pro b lem con sidered in this pa pe r ca n be form ulate d m athem atically as a c onstraine d optimization pro b lem: maximize PoV (x (I) ) (pos te rior value of inve n tory decision) s uch th at C x (I) ≤B (the inve n tory cos t is less tha n the bud get) where x (I) ∈ X (I) (the inve n tory methods a re chose n from ‘feasible’ s o lution s) , (3) with the notation exp l ained below. The “inve n tory de cision v e ctor” x (I) is a n in tege r v e ctor with n S ele me n ts. The s th ele me n t of x (I) , x (I) s ∈ { 1 , 2 , . . . , n I } , indicates the decision of the inve n tory method for tract s . We ass ume th at the inv e n tory methods a re “cos t-a ccura cy consis te n t”: Definition 1 Let I = { 1 , 2 , . . . , n I } be a set of i nv ent o ry m eth ods and let C s,i denote the cost of app l yi ng i nv ent o ry m eth od i to tract s . The set Iis cost-a cc ura cy co nsis t ent if: (1) C s,i > C s,j , i, j ∈ I, i  = j, imp l ies tha t fo r all s ∈ { 1 , 2 , . . . , n S } and for all a ∈ { 1 , 2 , . . . , n A } it holds σ2 a,s,i < σ2 a,s,j , (2) (witho ut lo ss of genera lity) the i nv ent o ry m eth o d in d ices a re ordered so that i > j, i, j ∈ { 1 , 2 , . . . , n I } imp l ies C s,i > C s,j for all s = { 1 , 2 , . . . , n S } . In other wor ds, mor e costly inve n tories a r e mor e accura te a nd the inve n tory methods h av e be e n orde re d s uch th at inv entory methods with hi ghe r index h av e hi ghe r a ccura cy. The feasib le s e t X (I) is formed by all inve n t ory decision vect ors in the s e t { 1 , 2 , . . . , n I } n S tha t sa tisfy the “pr ior var ia nce orde r” cons train t: Definition 2 An i nv ent o ry de cisio n ve ct o r x (I) = (x (I) 1 , . . . , x (I) n S ) follows prior va ri a nce order if for any two tracts s 1 , s 2 ∈ { 1 , 2 , . . . , n S } , s 1  = s 2 , it holds x (I) s 1 ≥x (I) s 2 when σ0 s 1 ≥σ0 s 2 , where σ0 s :=   n A a =1 (σ0 a,s ) 2 . This cons train t is a n assumpt ion of “rat ion ality” of the c osteffe ctiv e forest inv entories: for tw o inv e n tory methods a nd a ny tw o tracts, the inv e n tory tha t is mor e accura te can only be carried out for the tract with hi ghe r total prior unce rtain ty. By the indepe nde nce of the tracts and as s ortments for ( 1 ), the total preinve n tory s ta nda r d devia tion σ0 s of the timber volume can be c ompute d as in Definition 2. We wi l l d isc uss the pos te rior value PoV (x (I) ) in Section 5.1 . To incorporate the cos t-efficie ncy of the inve n t ory x (I) t o our formulation, w e m aximize PoV with respect to the fixed (monetary) inve n tory bud get B , which ca nnot be exc e e de d b y the fores t in- ve n tory cos t: C x (I) = n S  s =1 C s,x (I) s . (4) 5.1 The poste rio r value of the inve nto ry decisio n The pos te rior value PoV (x (I) ) measures on a n ave rage se nse the value of the decision x (I) , knowing that the inve n tory data Y | x (I) influences the future ha rves t scheduling decision. Thus, to define PoV (x (I) ) , we need to first define the future h arv est s cheduling pro b lem . The decision v ari ab le for the future h arv est sche duling deci sion i s an n T ×n S matrix X (T ) , with ele me n ts x (T ) ts ∈ { 0 , 1 } , t = 1 , 2 , . . . , n T , s = 1 , 2 , . . . , n S , de noting whethe r tract s is clear-cut at time t. The matrix X (T ) is s ubje ct to the cons train t  n T t=1 x (T ) ts ∈ { 0 , 1 } for all s = 1 , 2 , . . . , n S , since each tract is fully h arv este d exactly onc e or not at all. We de note b y X (T ) the s e t of all bin ary m atric es th a t sa ti sfy thi s cons train t. Furthe rmore, we use a dot ( ·) to denote “over all tracts”; the notation x (T ) t·, for example, refers to the row tof X (T ) . We model the efficiency of the future h arv est sche duling de cision in terms of a quadratic utility function: U(X (T ) , V) = − n T  t=1 n A  a =1  n S  s =1 x (T ) ts V a,s −D a,t  2 , (5) Downloaded from https://academic.oup.com/biometrics/article/80/3/ujae104/7774592 by Jyvaskyla University user on 27 September 2024 4 r Biometrics , 2024, Vol. 80, No. 3 which compares the total h arv est of as s ortment a at time tto its industry demand, D a,t , and penalizes devi ation s quadratically. The quadratic form not only facilit ate s the algorithms prese n ted in Section 6 but also aligns with the notion that the industry can adapt to small deviations from the de ma nd ta rgets but la rge devi ation s may caus e m a jor d iffic ulties . Note th at ( 5 ) is a random v ari ab le since it depends on the as s ortment vo lumes Vwhich are not known exactly before h arv esting. Conside r the n, how the inve n tory decision x (I) influe nces the future h arv est sche duling de cision X (T ) through the d ata o btained from the inve n tory decision: (1) When the inventory decision x (I) is made, a random in- ve n tory d atas e t Y ∼p Y| x (I) wi l l be o btained. The d ata Y are no is y measure me n ts of the unce rt ain volume s V. (2) Using the data Y, the kno wledg e reg arding Vcan be improv e d. This a moun ts to computing the pos te rior dis tri - bution of Vusing Bayes’ rule: p V| Y, x (I) (v | y, x (I) ) = p Y| V, x (I) ( y | v, x (I) ) p V ( v)  V p Y| V, x (I) ( y | v, x (I) ) p V ( v)d v . (3) Using the improv e d knowle dge of V, a better future harv est sche duling de cision X (T ) can be m ade. Ac c ording to Bayesian decision theory (cf. Raiffa and Schlaiffer, 1967 ; Hirshleife r a nd Riley, 1979 ), X (T ) should be chose n to ma ximize ex pected utility: max X (T) ∈X (T) E V| Y= y, x (I) [U(X (T ) , V)] = max X (T) ∈X (T)  V U( X (T ) , v) p V| Y, x (I) ( v | y, x (I) )d v. The qua n tity ca pturin g the logic of the abo ve steps (1)–(3) is the pos te rior value of the inve n tory decision: PoV (x (I) ) =  Y max X (T) ∈X (T)  V U (X (T ) , v) p V| Y, x (I) ×(v | y, x (I) )d v  p Y| x (I) (y | x (I) )d y, (6) which contains the outer integral over p Y| x (I) , since the data Y obt ained in ste p (1) a re ra ndom a nd unknown in the p l anning of the inve n tory decision x (I) . PoV (x (I) ) is r ela t ed t o the value of (imperfe ct) inform ation (VoI) (cf. Eidsvik et al., 2015 ) associated with x (I) : VoI (x (I) ) = PoV (x (I) ) −PV . (7) Here, the constant PV = max X (T) ∈X (T)  V U( X (T ) , v) p V ( v)d v(8) is the prior value, which corresponds to the m axim al expe cte d value gained from making the optimal future ha rves t scheduling de cision base d on p V (without c onside ring a ny inve n tory data c olle ction). The the oretical upper limit for ( 7 ) is the value of perfe ct inform ation obtaine d by rep l a cin g PoV (x (I) ) b y the pe rfect pos te rior v alue PoV ∗(s e e (2) in Web Appendix 7 ), th at is, the value of knowing volume Vexactly before the utility is maximized. 5.2 A co mputatio nall y att ractive rep res e ntatio n fo r the poste rio r value of inve nto ry decisio n Next we wi l l focus on a n alte rn ativ e re pre se n tation for ( 6 ) that allows us to desi gn a n al gorithm to s o lve the p l anning pro b lem ( 3 ) and to estimate VoI. To begin, note that the inne r in tegral in ( 6 ) can be simplified as follows:  V U( X (T ) , v) p V| Y, x (I) ( v | y, x (I) )d v = E { U(X (T ) , V) | Y = y, x (I) } = − n T  t=1 n A  a =1 E   n S  s =1 x (T ) ts V a,s −D a,t  2 Y a ·= y a ·, x (I)  = − n T  t=1 n A  a =1 ( x (T ) t·)  + a ·x (T ) t·+ ( x (T ) t·)  μ+ a ·−D a,t 2 , (9) where (X )  stands for the m atrix/v e ctor tran spos e of X , and w e h av e use d the formulas E [ X 2 ] = Var (X ) + E [ X ] 2 , and Var (a  Z ) = a  Cov (Z ) a for a random v e ctor Z . Furthermore, we use μ+ a ·:= E [ V a ·| Y a ·= y a ·, x (I) ] and + a ·:= Cov (V a ·| Y a ·= y a ·, x (I) ) , t o denot e the post erior mean and cov ari ance of V a ·, respe ctiv ely. Under the as sumption s of prior independence and conditionally indepe nde n t o bs erv ation s in the models of Section 4 , μ+ a ·is a v e ctor with ele me n ts E [ V a,s | Y a,s = y a,s , x (I) s ] , s = 1 , 2 , . . . , n S , and + a ·is a diagonal matrix with diagonal eleme n ts Va r (V a,s | Y a,s = y a,s , x (I) s ) , s = 1 , 2 , . . . , n S . To evaluate μ+ a ·and + a ·, it is ther efor e su ffic ient to compute only these s cal a r qua n t it ies, which can be done, for example, b y n ume rical integration. In Web Appendix 6 , it is shown that the maximization in ( 6 ) can be expres s ed as: min 1 2 n T  t=1 ( x (T ) t·)  Qx (T ) t·+ n T  t=1 ( c t )  x (T ) t·+ r, s.t. 0 ≤ n T  t=1 x (T ) ts ≤1 for all s = 1 , 2 , . . . , n S , where x (T ) ts ∈ { 0 , 1 } for all t, s, and Q : = 2 n A  a =1 + a ·+ μ+ a ·( μ+ a ·)  , c t := −2 n A  a =1 D a,t μ+ a ·, and r := n T  t=1 n A  a =1 D 2 a,t . (10) Thi s formulation allow s us t o writ e the post erior value ( 6 ) in the following form, which can be appr oxima ted b y Mon te Ca rlo Downloaded from https://academic.oup.com/biometrics/article/80/3/ujae104/7774592 by Jyvaskyla University user on 27 September 2024 Biometrics , 2024, Vol. 80, No. 3 r 5 s amp ling from the distribution p Y| x (I) : PoV (x (I) ) =  Y  −1 2 n T  t=1 ( ˆ x (T ) t·)  Q ˆ x (T ) t· − n T  t=1 (c t )  ˆ x (T ) t·−r  p Y| x (I) (y | x (I) )d y. (11) He re, ( ˆ x (T ) t·) 1 ≤t≤n T de notes the s o lution of ( 10 ), which depends on Q , the c t ’s and r(defined in ( 10 )). Note that Q and the c t ’s in turn depend on ythrough μ+ a ·and + a ·. 6 ALG ORITH MS FOR APPROXIMATING THE POSTERIOR VALUE AND SOLVING THE PLANNING PROBLEM The direct s o lution of the p l anning pro b lem ( 3 ) is d iffic ult, since the evaluation of the pos te rior value of inve n tory decision ( 11 ) r equir es the s o lut ion of a binary opt imizat ion pro b lem (the future h arv est sche duling problem) within a m ultiva riate in tegral. To ease this c omputation al pro b le m, we (1) a ppr oxima te the s o lut ion of the binary opt imizat ion pro b lem using a heuristic method and (2) appr oxima t e the int egral ( 11 ) by Monte Ca rlo sa mpling whe r e Mr ealiza tions ar e simula ted fr om dis tri - but ion p Y| x ( I) . Next, we wi l l prese n t thi s k ind of appr oxima tion for PoV (x (I) ) in Section 6.1 and then use it to construct an algorithm to s o lve the p l anning pro b lem ( 3 ) in Section 6.2 . 6.1 Appr oxim at ing th e poste rio r value The opt imizat ion pro b le m ( 10 ) within ( 11 ) ca n be see n as a speci al cas e of the G QA P, firs t proposed b y Lee a nd Ma ( 2004 ). Fol - lowing the formulation of the GQAP by Hahn et al. ( 2008 ), the pro b lem con siders the p l acement of Kfacilities to L locations with known quadratic a nd linea r coefficie n ts C ijkn a nd B ij , respe ctiv ely, and with constraints on the “space limit ations ”S j for each location j. In our pro b lem s e tting, the forest tracts are analogous to the facilities, and the locations to the times at which the tracts should be h arv este d. How ev er, in our s e tting, the G QA P is simplifie d s uch th a t ther e ar e no space limita tions, tha t is, S j = ∞ , and a large number of the coefficients B ij and C ijkn in the obje ctiv e function are zero. The G QA P i s NP-hard, as it i s a ge ne ralization of the ge ne ralized assi gnme n t pro b le m (Koopma ns a nd Be ckm ann, 1957 ) shown to be NP-hard by Fisher et al. ( 1986 ). Most of the research on the solution of the G QA P h as focuse d on exact solution using branch and bound methods (cf. Hahn et al., 2008 ; Pes s oa e t al., 2010 ). In our context, thes e k ind s of method s are too slow, since the dimension of the problem is r ela tively high, and w e ne e d to be able to s o lve the pro b lem Mtimes in order to evaluate one inve n tory decision, using the Monte Carlo approximation for ( 11 ). Th us, we rely ins tead on a fas t a nd gre e dy heurist ic method , which we call the “random sw e ep me thod”. The me thod con sid- e rs m ultiple ra ndom init ializat ions (s ta rting values) for the decision v ari ab le X (T ) . For each random init ializat ion, the method applies a se quenc e of “local moves” where each move improves the o bj e ctiv e function value. Onc e the s o lution cannot be improv e d further, the final o bj ect ive funct ion value is re c orde d and the pr ocess r epea ts for the next random init ializat ion. Finally, the me thod re turn s the s malle st o bj e ctiv e function value found among all random init ializat ions c onsidere d. In our context, this strategy is computationally a t tractive for two reas on s. First, the change in the o bj ect ive funct ion value for each local move is cheap to evaluate in comparison to the full o bj e ctiv e function value of ( 10 ). Se c ond, the local mov es can explore the full feasible region and always preserve the feasibility of the impro vin g candid ate s o lution s. The local moves correspond to changes in the ha rves t schedul - ing of a single tract while the h arv est sche duling of the other tracts remains fixed. Invokin g ea ch local move changes the obje ctiv e function value of ( 10 ) by an “adjustment”. For each tract s , the re a re n T + 1 pos sib le adjus tme n ts: the adjus tme n t corresponding to the local move that does not change the timing (0) and the adjustments for mov es th at alter the h arv est timing (give n curre n t timing): From time u to ‘no h arv est’: −c us −(Q ·s )  x (T ) u ·+ 1 2 Q ss , (12) From ‘no h arv est’ to time u : c us + (Q ·s )  x (T ) u ·+ 1 2 Q ss , (13) From time u to time v : −c us + c vs + (Q ·s )  x (T ) v · −(Q ·s )  x (T ) u ·+ Q ss , (14) whe re c us s ta nds for the s th ele me n t of the v e ctor c u , Q ss is the s th diagonal ele me n t of Q , a nd Q ·s is the s th column of Q . The adjus tme n ts in ( 12 )–( 14 ) h av e be e n de riv e d by inspe cting the chan gin g terms in the obje ctiv e function when the local moves a re a pplied. Usin g a djus tme n ts ( 12 )–( 14 ), Al gorithm 1 give s pse udo-code for the random sw e ep method th at can be use d to a pproxi - mately s o lve the opt imizat ion pro b lem ( 10 ) given the spec i fication ( Q , ( c t ) 1 ≤t≤n T , r) and n ri random init ializat ions. The name RandomSweep for the algorithm comes from line 7, where a random processing order (a “sw e ep” ov er the tracts) is drawn. Lines 6–15 carry out the core computa tion, wher e local mo ves and a djustments are chosen for the tracts, and the s o lution and the current o bj e ctiv e function value are updated , unt il the curre n t s o lution cannot be improv e d further. On line 10, there may be multiple local moves with ne gativ e adjus tme n ts, a nd the choice of the local move in such a situation affects how greedy the resulting algorithm wi l l be. The greediest choice is to always select the local move that improves (decreases) the obje ctiv e function value the most. In this pa pe r, we refer to this strategy as the “gre e dy desc e n t”s trategy. Ins tead, it is als o pos sib le to choos e a les s gre e dy desc e n t s trate gy th at w e call “con s erv ative des ce n t ”, w hich always chooses the local move that improves the o bj ect ive funct ion value the least (but improves nonethe le ss). The c onservativ e desc e n t s tra tegy incr eases the Downloaded from https://academic.oup.com/biometrics/article/80/3/ujae104/7774592 by Jyvaskyla University user on 27 September 2024 6 r Biometrics , 2024, Vol. 80, No. 3 Al go r ithm 1 RandomSweep ( Q , ( c t ) t=1 , 2 , ... ,n T , r;n ri ) 1: Initi alis e the return value m = ∞ . 2: for i ∈ 1 , 2 , . . . , n ri do 3: Sample a random feasible decision matrix X cur . 4: Set m cur to the full obje ctiv e function value associated with X cur . 5: If m cur < m , s e t m = m cur , X best = X cur . 6: while true do 7: Sa mple a ra ndom pr ocessing or de r p (a pe rm ut ation of tract indice s 1 , 2 , . . . , n S ). 8: for tract index s ∈ p do 9: With X cur and (Q , c t ) , compute all adjustments for tract s using (12)-(14). 10: Choose a local move with adjus tme n t c ≤0 . De note b y l the chose n local move. 11: Apply local move l to the s th column of X cur . 12: Set m cur = m cur + c . If m cur < m , s e t m = m cur , X best = X cur . 13: en d fo r 14: If X cur did not change in the above loop, break out of the while loop. 15: end while 16: en d fo r 17: Return m , X best . c onv e rge nce time of the method. How ev er, ev en with this strat- e gy, w e h av e o bs erv e d th at the while loop in Algorithm 1 very ofte n s ti l l terminates in at most a few dozen iterations. We wi l l inves ti gate the effects of the choice of the desce n t s trategy in Section 6.3 . We re ma rk that Al gorithm 1 is similar to the Hero method of P ukkala a nd Ka ngas ( 1993 ), which h as be e n found to s trike a n efficie n t tradeoff in terms of c omputation al c omplexity vs . quality of s o lut ion (Puk kala and Kurtt il a, 2005 ) in forest p l anning pro b lem s. Algorithm 1 extends Hero by introducing the random pr ocessing or der of the tracts on each iteration; the Hero me thod would in stead proces s the tracts without v arying the orde r. This ra ndomization e ns ures th at no pa rticula r tract has “priority” over the other tracts and introduces v ari ation to the s o lution s exp lored by the me thod. Furthe rmore, in con tras t to He ro, Algorithm 1 provides “gre e dyness c ontrol” with the tw o desc e n t s trategies a nd has bee n tailored to s o lve ge ne ral opt imizat ion pro b lem s of the form ( 10 ) by taking adva n tage of the adjus tme n t formulas ( 12 )–( 14 ). Note that it is not strictly ne c es s ary to use Algorithm 1 to s o lv e the h arv est sche duling proble ms a rising in this work; we d isc uss refine me n ts a nd alte n ativ es in Se ction 8 . Al go r ithm 2  PoV ( x (I) ;M, n ri ) 1: for k = 1 , 2 , . . . , M do 2: Simula te y (k) fr om p Y| x (I) . 3: Given y (k) , compute μ+ a ·and + a ·for each a = 1 , 2 , . . . , n A . 4: Using ( μ+ a ·, + a ·) 1 ≤a ≤n A a nd (D a,t ) 1 ≤a ≤n a , 1 ≤t≤n T , cons truct Q , (c t ) 1 ≤t≤n T a nd r in (10). 5: m (k) , X (k) best ← RandomSweep (Q , (c t ) 1 ≤t≤n T , r;n ri ) 6: en d fo r 7: return −1 M  M k=1 m (k) With the help of Algorithm 1, we prese n t Al gorithm 2 which is a Mon te Ca rlo a ppr oxima tion of ( 11 ) give n the (ca ndidate) inve n tory decision x (I) , MMonte Carlo s amp les and n ri random init ializat ions for Algorithm 1. For each ite ration k, Al gorithm 2 proc e e ds b y sim ula ting a measur e me n t v e ctor y (k) from the m argin al distribution p Y| x (I) (line 2). This can be carried out using the models of Section 4 by first simulating vfrom the prior dis tribution p V , a nd the n sim ulat ing y (k) ∼p Y| V, x ( I) (·| V = v, x (I) ) . The measure me n t y (k) is s ubse que n tly use d to c ompute the pos te rior mea ns a nd cova ria nces ( μ+ a ·, + a ·) 1 ≤a ≤n A (line 3). Finally, these qua n t it ies (to ge ther with the de ma nds D a,t ) a re use d to c onstruct the minimiza tion pr oblem ( 10 ), which is then appr oxima tely s o lved using Algorithm 1 and the minimum m (k) re c orde d (line 5). The use of Algorithm 1 on line 5 can be rep l ac e d with a nothe r opt imizat ion me thod targe ting the pro blem ( 10 ). The output of Algorithm 2 is the mean of the ne gativ e m (k) ’s, where the change of sign is applied since we are int erest ed in maximizing rather than minimizing. Fin ally, w e rem ark th at with minor modifications to the minimiza tion pr o b le m ( 10 ) a nd Al gorithms 1 a nd 2, the methods in this se ction m ay also be lev erage d to approxim ate PV and PoV ∗. The details are given in Web Appendix 7 . 6.2 Th e pract ic a l s olutio n m eth od Algorithm 3 pres ents ps eudo-code for our practical method that appr oxima tely s o lves the p l anning pro b le m ( 3 ) for a n inve n tory budget B , using MMonte Carlo s amp les and n ri random initializations. Al go rithm 3 So lv ePlannin gPro blem (B ;M, n ri ) 1: Initi alis e an empty container O for res ults . 2: Construct the set X (I) . 3: for x ∈ X (I) do 4: Evaluate p =  PoV (x ;M, n ri ) (Algorithm 2). 5: Evaluate c = C x using (4). 6: Save (x, p, c ) to the con taine r O . 7: en d fo r 8: Set O B = { (x, p) : (x, p, c ) ∈ O and c ≤B } 9: Find (x ∗, p ∗) from O B such that p ∗≥p for all ele me n ts (x, p) ∈ O B . 10: return x ∗, p ∗ The al gorithm firs t cons tructs the feasib le s e t X (I) , a nd the n, for each x ∈ X (I) , computes the inve n tory cos t using ( 4 ) and the appr oxima t e post erior value using Algorithm 2. The output of Algorithm 3 is the highest appr oxima t e post erior value and the as s oci ate d inv e n tory decision with cost less than B (lines 8– 9). The inve n tory bud get B is only r equir ed by Algorithm 3 after processing the set X (I) . Ther efor e, Algorithm 3 can s o lve ( 3 ) appr oxima tely for multip le budge ts with a single run of the processing loop (lines 3–7). Under the assumption of “cost-a ccura cy consis te ncy” (Defini - tion 1) and distinct total prior volume va ria nc es, each inv e n tory de cision v e ctor x ∈ X (I) can be c onstructe d as follows . First, define κ(n I , n S ) := { (z 1 , z 2 , . . . , z n S ) : z i ∈ { 1 , 2 , . . . , n I } , z i ≤z j ∀ i < j} , that is, κ(n I , n S ) is the s e t of or der e d v e ctors of length n S with ele me n ts in { 1 , 2 , . . . , n I } . For each z ∈ κ(n I , n S ) , we can construct a corresponding x ∈ X (I) by s e tting x o j = z j for j ∈ 1 , 2 , . . . , n S , where ois the v e ctor of tract indic es sorte d in the order of increasing total prior s ta nda r d devia tion σ0 s (see Definition 2). It is well known that | κ(n I , n S ) | = n S + n I −1 n S , implying |X (I) | = n S + n I −1 n S , as well. 6.3 Inv estig ating the pe rfo rma nce of Alg o rithm 1 Our appr oxima tion of the pos te rior v alue us es Algorithm 1 to s o lve the minimization pro b le m ( 10 ). The pe rforma nce of the al gorithm is the r efor e importa n t for the ove rall pe rforma nce of Downloaded from https://academic.oup.com/biometrics/article/80/3/ujae104/7774592 by Jyvaskyla University user on 27 September 2024 Biometrics , 2024, Vol. 80, No. 3 r 7 FIGURE 1 The variability in the final obje ctiv e function values of Algorithm 1 with respect to the n umbe r of random init ializat ions and the desce n t s trategy d isc ussed in Section 6.1 . For both desce n t s trategies, the depicted lines correspond to the 1% and 99% quantiles of final o bj e ctiv e function values obtained, when Algorithm 1 was run 10 000 times with the number of init ializat ions shown in the horizontal axis. Note that the horizontal axis does not increase linearly. The panes indicate the inventory decision for all 100 tracts that was used to construct the three minimization pro b lem s of the form ( 10 ) studied here. the appr oxima tion a nd th us the s o lution of the full p l anning pro b lem ( 3 ). As Algorithm 1 fea tur e s randomne ss, w e studie d the v ari ation in the found optim a. We c onstructe d thre e te st minimiza - tion pro b lem s of the form ( 10 ), as fo llows. First, usin g the tra ct data of Section 3 , we simulated three inve n tory d atas e ts with 5, 10, or 20 plots for every tract from the normal model of Web Appendix 5 . For each test pro b lem, we ran Algorithm 1 10 000 times using both desce n t s tra tegies fr om Section 6.1 ) a nd va ried the n umbe r of ra ndom init ializat ions n ri from 5 up to 20 000. After each run, we r ecor ded the final o bj ect ive funct ion value. Figure 1 s umm arize s the re sults of this expe rime n t b y disp l aying the empirical 1% and 99% qua n tiles of the final 10 000 o bj ect ive funct ion values given the desce n t s trategy a nd the n umbe r of random init ializat ions n ri of Algorithm 1. As expe cte d, the va ri - ability of the final o bj e ctiv e function values quickly diminishes as n ri increases. As the opt imizat ion pro b le m is NP-ha rd, it is pos sib le that the true global minimum is even lower than the achieved s o lution s. To study whe ther there exist better s o lution s, we further ran the Gurobi optimiser (Gurobi Opt imizat ion, LLC, 2023 ) for 14 h in the case of the minimiza tion pr o b lem con structed when the inventory decision was 20 plots for each tract. Our imple me n tation in Gurobi used Algorithm 1 as a heuristic for the mixed -in tege r s o lver. The results for this expe rime n t a re shown in Web Figure 4 in Web Appendix 8 and indicate that the minima found by Algorithm 1 are typically within 5%– 15% of the best minimum found using Gurobi when a sufficie n t n umbe r of ra ndom init ializat ions are used . Note that compared to the 14 h runtime of Gurobi, Algorithm 1 finishes in seconds. In te res tingly, in Fi gure 1 , the con s erv ative des ce n t s trategy exhibits les s v ari ability for all n ri s tudied, a nd also a ppea rs to find s o lution s with sli gh tly bette r o bj e ctiv e function values . Our intuition is that this occ ur s since the con s erv ative des ce n t exp lores the s o lution space more efficie n tly b y avoidin g mo ves tha t impr ove the o bj ect ive funct ion value too e ag e rly, whe n the ra ndomly sa mpled solution candidates wi l l sti l l be of poor quality. Finally, the panes f rom lef t to right in Figure 1 show decreasing o bj e ctiv e function values as the inve n tory deci sion i s made more accurate and expensive. This makes s en s e, since more accurate inve n tory data wi l l provide better information on the unknown volumes V, impro vin g the found o bj e ctiv e function value in the h arv est sche duling pro b lem ( 10 ). Downloaded from https://academic.oup.com/biometrics/article/80/3/ujae104/7774592 by Jyvaskyla University user on 27 September 2024 8 r Biometrics , 2024, Vol. 80, No. 3 FIGURE 2 The inve n tory decisions returned by Algorithm 3 for the Sw e dish tract d atas e t with budget cons train ts from 10 000 to 25 000 by incre me n ts of 1000. The shade of a cell shows the pr eferr ed inve n tory decision for a tract (vertical axis) under the given budget (horizontal axis). 7 APPLICATION TO THE SWEDISH TRACT DATASET Next, we apply the methods of Section 6 to the Sw e dish tract d atas e t dis cus s e d in Se ction 3 . For e ach budg e t, we app ly Algorithm 3 with the model d isc usse d in Se ction 4 . Thi s yield s the pr eferr e d inv e n tory me thod (s amp ling 5, 10, or 20 p lots) for ea ch tra ct in the d atas e t. Motiv ated by our findings in Section 6.3 , we use 2500 Mon te Ca rlo sa mples, 2500 ra ndom init ializat ions and the con s erv ative des ce n t s trategy for Al gorithm 1 within Al - gorithm 3. Figure 2 s umm arizes the res ults . The inv e n tory decision s fo llow an organized structure due to the ordering of the tracts by the prior volume va ria nce (on the vertical axis ). F or the extreme bud gets 10 000 a nd 25 000, the chea pes t (5 plots inve n tory for every tract) and the most expensive inve n tory (20 plots inventory for every tract) are chosen, as the costs of these inventories coincide exactly with the extreme budgets . In c ontrast, all thre e inv e n tory methods a re us ed for s ome tracts when the budge t is be tw e e n 11 000 a nd 23 000. For these budgets, the proportions of the inve n tory methods sele cte d for each tract vary, which occ ur s because of the diffe re nce betwee n the as s ortme n t composition in the tracts and the diffe re n t cos ts of the inve n tories. For budgets 16 000–19 000, it a ppea rs tha t the pr oportions of tracts with 10 or 20 plots inve n tory a r e r oughly the sa me, whe reas for budgets 20 000–25 000, the 20 plots’ inve n t ories begin t o “domina te”, as ther e is enough money to spend in the bud get. For bud - ge ts les s th an 15 000, espe ci ally, the proportion s of tracts to measure with 10 and 20 plots vary without a clear pa t tern. Figure 3 disp l ays the m axim al pos te rior v alue as s oci ated with the decisions computed by Algorithm 3 th at c oste d less than each inve n tory bud ge t con sidered . In addit ion, the prior value PV and the perfect posterior value PoV ∗d isc ussed in Section 5.1 are shown with horizontal lines. The increasing pos te rior values depict the increasing information gain from la rge r bud gets. In comparison to PV , the figure shows that a si gnifica n t information gain is obtaine d ev e n b y ca rrying out the chea pes t pos sib le fores t inve n tory cos ting 10 000. The dis crepancy be tween PoV ∗ and the m axim al pos te rior v alue o btained for budge t 25 000 depicts the information gain that could sti l l be obtained by even more expensive and accurate inve n tory methods. Downloaded from https://academic.oup.com/biometrics/article/80/3/ujae104/7774592 by Jyvaskyla University user on 27 September 2024