A Node Formulation for Multistage Stochastic Programs with Endogenous Uncertainty
AA Node Formulation for Multistage Stochastic Programs withEndogenous Uncertainty
Giovanni PantusoDepartment of Mathematical SciencesUniversity of Copenhagen, 2100 Copenhagen, [email protected] 16, 2021
Abstract
This paper introduces a node formulation for multistage stochastic programs with endogenous (i.e.,decision-dependent) uncertainty. Problems with such structure arise when the choices of the decisionmaker determine a change in the likelihood of future random events. The node formulation avoidsan explicit statement of non-anticipativity constraints, and as such keeps the dimension of the modelsizeable. An exact solution algorithm for a special case is introduced and tested on a case study. Resultsshow that the algorithm outperforms a commercial solver as the size of the instances increases.
Multistage stochastic programs offer a viable framework for modeling and solving problems involving asequence of decisions interspersed with partial resolutions of some stochastic process. At each decisionstage the decision maker knows the content of the uncertainty resolved until that stage and a probabilisticcharacterization of the remaining stochastic process. In the classical settings, it is assumed that decisions donot modify the stochastic process in any way (see e.g., Kall and Wallace, 1994; Birge and Louveaux, 1997).In other words, the uncertainty is entirely exogenous. While this description fits a large number of decisionproblems, several other can be found where decisions have an influence on the remainder stochastic process,e.g., by changing the likelihood of future realizations. This category of problem is referred to as multistagestochastic programs with endogenous uncertainty (MSPEU).Following (Goel and Grossmann, 2006), there exist at least two ways in which decisions can influencethe underlying stochastic process. The first possibility is that decisions alter the probability distribution ofthe stochastic process, thus changing the likelihood of the possible events. The second possibility is thatdecisions determine the time when the uncertainty is (partially) resolved. This article is concerned withmultistage stochastic programs affected by the first type of endogenous uncertainty.The research dealing with decisions influencing probability distributions is rather sparse. Jonsbr˚atenet al. (1998) consider a case where decisions influence both the probability measure and the timing ofthe observation, i.e., at which stage certain random variables will be observed. The framework includesboth two-stage and multistage problems, though the decisions influencing the uncertainty must be madeat the first-stage only. Ahmed (2000) illustrates examples of problems with endogenous uncertainty, suchas facility location, network design and server selection. The author presents an exact solution methodfor the resulting one-stage integer problems. Viswanath et al. (2004) consider the problem of investing instrengthening actions for the links of a network subject to disruptive events. The problem is modeled asa two-stage stochastic program where first-stage investment decisions influence the likelihood of disruptiveevents happening at upgraded links. The same problem is studied also by da Costa Flach (2010), Peetaet al. (2010) and Laumanns et al. (2014). Held and Woodruff (2005) consider the problem of interdictinga stochastic network, that is a network whose structure is unknown to the interdictor. In this problem,the probabilities of different future network configurations depend on previous interdiction actions. Tonget al. (2012) present an oil refinery planning problem considering that the uncertainty in product yield is1 a r X i v : . [ m a t h . O C ] F e b nfluenced by operation mode changeovers. Hellemo (2016) and Hellemo et al. (2018) discuss several waysof incorporating the influence of decision variables on the underlying probability distributions in two-stagestochastic programs. Particularly, the authors formulate two-stage models where prior probabilities aredistorted through an affine transformation, or combined using a convex combination of several probabilitydistributions. Furthermore, the authors present models which incorporate the parameters of the probabilitydistribution as first-stage decision variables. Finally, Escudero et al. (2018) study the problem of mitigatingthe effects of natural disasters through preventive actions. The problem is formulated as a three-stagestochastic bilinear integer program with both exogenous and endogenous uncertainty. Particularly, decisionscan influence both the probabilities and the intensity of future uncertain events.When the second type of uncertainty is considered, the decision making process is such that the un-certainty is not resolved automatically at each decision stage as in classical stochastic programs with onlyexogenous uncertainty. Rather, the decisions made implicitly determine the time when the uncertainty isresolved. A typical example, provided by Goel and Grossmann (2004) and based on the petrochemical in-dustry is as follows. A decision maker has to decide which gas reservoirs to explore, and when, by installingexploration facilities. The size and quality of the reservoirs is uncertain and can be known only when fa-cilities have been installed. Thus, the time the uncertainty is resolved depends on installation decisions.The literature dealing with this type of stochastic programs includes Colvin and Maravelias (2008), Tarhanand Grossmann (2008), Tarhan et al. (2009), Colvin and Maravelias (2010), Gupta and Grossmann (2011),Mercier and Van Hentenryck (2011), Tarhan et al. (2013), Apap and Grossmann (2016).The contribution of this paper is as follows. First, we introduce a novel node-formulation for multi-stage stochastic programs where decisions (at all stages) influence probability distributions at future stages.A node formulation automatically ensures non-anticipative decisions and thus avoids writing explicit non-anticipativity constraints (NACs) that are necessary when the uncertainty is represented via scenarios. NACshave typically been addressed as a bottleneck of available models for MSPEU, and motivated recent researchto find and discard redundant NACs, see Apap and Grossmann (2016), Hooshmand and MirHassani (2016)and Hooshmand and MirHassani (2018). The new formulation is based on a novel scenario tree structurewhich incorporates the possibility that several (finitely many) different distributions for later stages can em-anate based on the decisions made at a certain decision stage. The new scenario tree structure represents thesecond contribution of this paper. Third, we propose a general-purpose efficient algorithm for a special classof MSPEUs. We demonstrate the use of the new formulation, scenario tree structure and solution algorithmon instances of the Football Team Composition Problem . The instances are made available online for thebenefit of future research at https://github.com/GioPan/instancesFTCPwithEndogenousUncertainty .The remainder of this paper is organized as follows. In Section 2 we provide a model formulation andnew scenario tree structure for MSPEUs. In Section 3 we describe a solution algorithm for a special caseof MSPEUs. In Section 4 we present a case study in which we formulate the Football Team CompositionProblem as a MSPEU and solve it using our specialized algorithm. Finally, we provide concluding remarksin Section 5.
The decision maker is concerned with a sequence of decisions ( x t ) Tt =1 at decision stages t = 1 , . . . , T , condi-tional on the partial resolution of a random process ( ξ t ) Tt =1 . At decision stage t , decisions are non-anticipative,meaning that they are based only on the realization of the random process up to, and including, ξ t . Therealization of the remaining random process ξ t +1 , . . . , ξ T is still uncertain. For t = 2 , . . . , T , the probabilitydistribution of ξ t , . . . , ξ T is dependent on past decisions x , . . . , x t − . The resulting multistage stochasticprogram is thus characterized by endogenously defined uncertainty.Consistently with Jonsbr˚aten et al. (1998) we assume that the set of potential probability distributionsenforced by decisions is finite and countable. Furthermore, we assume that the probability distributionsare discrete (possibly after a scenario generation phase). The latter assumption is rather standard andrequired in order to solve real-life stochastic programs (except perhaps for a number of specific applications).The resulting decision-dependent discrete stochastic process can be depicted by means of the scenario treestructure illustrated in Figure 1. In what follows we refer to this scenario tree structure as a multi-distribution T N ,d d ∈ D N md d ∈ D m m N m,d +1 d + 1 ∈ D m lqo N nd d ∈ D n n N n,d +1 d + 1 ∈ D n T − τ N ,d +1 d + 1 ∈ D Figure 1: Multidistribution scenario tree. scenario tree (MDST) (see, e.g., also how Kaut et al. (2014) modified the classical structure of scenario treesin order to account for multiple time resolutions).In an MDST the root node, arbitrarily named 0, represents the current state of the world, when first-stagedecisions are made. First-stage decisions will enforce one out of a finite number of distributions representedby the set D . As an example, in Figure 1, decisions might determine, among other, distributions d ∈ D or d + 1 ∈ D . Distribution d is characterized by realizations, represented by nodes, that include l and m , whiledistribution d + 1 is characterized by realizations that include q and o . Similarly, at stage t = 2, assumingrealization m of distribution d occurs at stage t = 1, the actions of the decision maker will determine oneout of a number of different distributions including d and d + 1 from the set D m . This process continues ina similar manner until stage T −
1. Given a realization n at stage T − D n , which include d and d + 1. Finally, at stage T all the uncertainty is resolved and the decision maker makes final decisions.Given an MDST, let us introduce the notation necessary for formulating the MSPEU. Let N be the setof nodes in the scenario tree, N t be the set of nodes at stage t , 0 the root node, t ( n ) ∈ { , . . . , T } the stageof node n , and a ( n ) ∈ N the parent node of node n except the root node. Let D n be the set of possible3 N ,d d ∈ D N ,d d ∈ D N ,d d ∈ D N ,d d ∈ D N ,d d ∈ D N ,d d ∈ D N ,d d ∈ D N ,d d ∈ D N ,d d ∈ D N ,d d ∈ D N = { , . . . , }D = D = . . . = D = { d , d }D = . . . = D = ∅N ,d = { , }N ,d = { , }N ,d = { , }N ,d = { , }N ,d = { , }N ,d = { , }N ,d = { , }N ,d = { , }N ,d = { , }N ,d = { , } Figure 2: Example multidistribution scenario tree with three stages, two possible distributions for each node,and two possible realizations for each distribution.distributions which can be enforced by decision made at node n ∈ N , and N nd be the child nodes of node n if distribution d ∈ D n is enforced. Let π n be the probability of node n with π = 1, and (cid:80) m ∈N dn π m = π n for all n ∈ N , d ∈ D n . An example of this notation is provided in Figure 2 for a three-stage scenario tree,with two possible distributions emanating from each node and each distribution being characterized by twopossible realizations.Let decision variables x n ∈ R N t ( n ) , n ∈ N represent decisions made at node n . These decisions maybe integer or fractional and represent ordinary decision made in the decision process. Let δ nd be a binaryvariable which captures the probability distribution enforced by the decisions made at node n . It takes value1 if the decisions at node n enforce probability distribution d at the child nodes of n , 0 otherwise. Finally, let θ n ∈ R be a decision variable which holds the expected value of the decisions made at the nodes descendingfrom n . For the reader’s convenience, the notation is also reported in Appendix A in a tabular format. AnMSPEU is thusmax r T x + (cid:88) d ∈D q d δ d + θ (1a)s.t. (cid:88) d ∈D n δ nd = 1 n ∈ N , (1b)4 n x n + (cid:88) d ∈D n B nd δ nd + C a ( n ) x a ( n ) + (cid:88) d ∈D a ( n ) D a ( n ) ,d δ a ( n ) ,d = h n n ∈ N , (1c) θ n = (cid:88) d ∈D n δ nd (cid:18) (cid:88) m ∈N nd π m ( r Tm x m + (cid:88) d ∈D m q md δ md + θ m ) (cid:19) n ∈ N \ N T , (1d) θ n = Θ n n ∈ N T , (1e) x n ∈ X t ( n ) n ∈ N , (1f) δ nd ∈ { , } n ∈ N , d ∈ D n , (1g) θ n ∈ R n ∈ N . (1h)where, for each n ∈ N , r n ∈ R N t ( n ) and q nd ∈ R represent the rewards of decisions x n and δ nd ,respectively, A n ∈ R M t ( n ) × N t ( n ) , B nd ∈ R M t ( n ) × , C a ( n ) ∈ R M t ( n ) × N t ( a ( n )) , and D a ( n ) ,d ∈ R M t ( n ) × arematrices of coefficients and h n ∈ R M t ( n ) a right-hand-side vector, with the assumption that C a (0) ,d := and D a (0) ,d := . Finally, Θ n ∈ R represents the future expected value at leaf node n , and X t ( n ) ⊆ R N t ( n ) thedomain of the x n variables. Objective function (1a) represents the sum of the profit for the decisions made athe root node ( n = 0) and expected profit of future decisions. Constraints (1b) ensure that the decisions madeat each node determine exactly one probability distribution. That is, the stochastic phenomena followingthe decisions at node n , are described by exactly one probability distribution (among the available ones),enforced by the decisions made. Constraints (1c) describe the dependency between decision variables x and δ at each node, and between them and the corresponding decision variables at the parent node. That is,the choice of a probability distribution at node n , δ nd , depends on the decisions x n made at the same nodeas well as on the decisions x a ( n ) made at the parent node and on the consequent choice of a probabilitydistribution, δ a ( n ) ,d . Constraints (1d) ensure that decision variables θ n hold the expected value of futuredecisions calculated according to the probability chosen. Consider a generic node n other than a leaf node,and note that, according to (1b), at this node there will be exactly one δ nd equal to one, that is exactlyone distribution be chosen. Thus, the right-hand-side of constraints (1d) will be equal to the term of theouter summation corresponding to the index d whose δ nd is set to one. The remaining terms are zero.Correspondingly, θ n will hold the expected value calculated according to the chosen probability distribution d . Constraints (1d) can be linearized using standard techniques. A linear reformulation is provided inAppendix B and a general procedure to determine the necessary big- M values is described in Appendix C.Constraints (1e) set the future expectation at the leaf nodes. Finally, constraints (1f) to (1h) set the domainof the decision variables. Particularly, X t ( n ) represents the domain of the x n variables and may imposeintegrality restrictions on some/all variables.A node formulation implicitly includes non-anticipativity, that is, automatically ensures that the decisionsmade at a given stage are only based on available information. On the other hand, a classical scenarioformulation requires non-anticipativity explicitly enforced by means of constraints, and in general, generatesa much larger problem. This can be illustrated by the following example. Consider Figure 3 which providesthe scenario representation of the example MDST in Figure 2. The number of scenarios is 15, and is thesame as the number of leaf nodes in Figure 2. Assuming the decision at each stage are represented by N decision variables and have to satisfy M constraints, the corresponding scenario formulation would include • N × × • M × × • approximately N × (15 + 8) non-anticipativity constraints (approximately 15 for the first stage and 8for the second stage, though more efficient specifications may be possible).The node formulation would include • N ×
21 variables (where 20 is the number of nodes in the scenario tree in Figure 2), and • M ×
21 constraints.Clearly, node formulations generate, in general, much smaller problems with non-anticipativity constraintsplaying an important role in scaling up the dimension of a scenario formulation, see e.g., Apap and Grossmann(2016), Hooshmand and MirHassani (2016) and Hooshmand and MirHassani (2018) for how to reduce thenumber of non-anticipativity constraints. 5
Figure 3: Scenarios in the example multidistribution scenario tree in Figure 2. Plain lines connect the nodesbelonging to the same scenario. Dashed lines represent non-anticipativity constraints.6
A solution algorithm for a special case
In this section we introduce an algorithm for solving the special case of problem (1) with C a ( n ) = for all n . The algorithm builds on the fact that, when C a ( n ) = , the link between decisions at subsequent stagesis created by the δ variables, i.e., those that enforce a probability distribution for the next stage. In thiscase, since we have finitely many distributions, we are allowed to enumerate the expected values obtainableat each node in the scenario tree. This procedure would become impractical when C a ( n ) is different from ,unless further assumptions on x n are made.Let us thus consider the following equivalent formulation of problem (1). z = max r T x + (cid:88) d ∈D q d δ d + θ (2a)s.t. (cid:88) d ∈D n δ nd = 1 n ∈ N , (2b) A n x n + (cid:88) d ∈D n B nd δ nd + (cid:88) d ∈D a ( n ) D a ( n ) ,d δ a ( n ) ,d = h n n ∈ N , (2c) θ n ≤ φ nd + M nd (1 − δ nd ) n ∈ N \ N T , d ∈ D n (2d) φ nd = (cid:88) m ∈N nd π m ( r Tm x m + (cid:88) d (cid:48) ∈D m q md (cid:48) δ md + θ m ) n ∈ N \ N T , d ∈ D n (2e) θ n = Θ n n ∈ N T , (2f) x n ∈ X t ( n ) n ∈ N , (2g) δ nd ∈ { , } n ∈ N , d ∈ D n , (2h) θ n ∈ R n ∈ N , (2i) φ nd ∈ R n ∈ N \ N T , d ∈ D n . (2j)Problem (2) modifies problem (1) in two elements. First, constraints (2c) take into account that C a ( n ) = .Second, constraints (1d) have been linearized using constants M nd and auxiliary decision variables φ nd yielding constraints (2d) and (2e) (see e.g., Appendix C for a general purpose procedure to determine thesebig- M values).Problem (2) can be solved in a dynamic programming fashion using the following backward procedure.The procedure starts by calculating the optimal last-stage expectation for each node at the second-last stageand for each possible probability distribution. For each m ∈ N T − and probability distribution k ∈ D m , theoptimal expectation at the last stage is obtained by solving the following problem.Φ mk = max (cid:88) n ∈N mk π n ( r Tn x n + (cid:88) d ∈D n q nd δ nd + Θ n ) (3a)s.t. (cid:88) d ∈D n δ nd = 1 , n ∈ N mk , (3b) A n x n + (cid:88) d ∈D n B nd δ nd = h n − D mk n ∈ N mk , (3c) x n ∈ X T n ∈ N mk , (3d) δ nd ∈ { , } n ∈ N mk , d ∈ D n . (3e)Problem (3) provides the last-stage expectation for each node m ∈ N T − at the second-last stage, and foreach possible selection of a probability distribution from D m . Notice in (3) that Θ n is input data and that D mk is moved to the right-hand-side to stress that it does not multiply a decision variable as in (2c), sincethe choice of a probability distribution at the parent node m has been fixed to k .Then, for stage t = T − , . . . ,
1, for node m ∈ N t , and for distribution k ∈ D m we calculate Φ km bysolving problem (4).Φ mk = max (cid:88) n ∈N mk π n ( r Tn x n + (cid:88) d ∈D n q nd δ nd + θ n ) (4a)7.t. (cid:88) d ∈D n δ nd = 1 , n ∈ N mk , (4b) A n x n + (cid:88) d ∈D n B nd δ nd = h n − D mk n ∈ N mk , (4c) θ n ≤ Φ nd + M nd (1 − δ nd ) , d ∈ D n , n ∈ N mk , (4d) x n ∈ X t ( n ) n ∈ N mk , (4e) δ nd ∈ { , } n ∈ N mk , d ∈ D n , (4f) θ n ∈ R n ∈ N mk . (4g)Notice in problem (4) that Φ nd is input data and has been calculated in the previous steps of the algorithm.Finally, we can solve the following problem for the root node. z = max r T x + (cid:88) d ∈D q d δ d + θ (5a)s.t. (cid:88) d ∈D δ d = 1 , (5b) A x + (cid:88) d ∈D B d δ d = h , (5c) θ ≤ Φ d + M d (1 − δ d ) , d ∈ D , (5d) x ∈ X , (5e) δ d ∈ { , } d ∈ D , (5f) θ ∈ R . (5g)Letting D = max n ∈N |D n | , the algorithm entails solving O (cid:0) D |N | (cid:1) mixed-integer programs of size signif-icantly smaller than (1). In this section we present a case study based on the
Football Team Composition Problem (FTCP – Pantuso(2017); Pantuso and Hvattum (2020)) which we extend in order to account for decision-dependent uncertainty.The problem consists of selecting players for a football team while their future market value is uncertainand influenced by the team for which they play. The scope of the computational study is to compare theperformance of the algorithm to that of a state-of-the-art commercial solver. We remark, however, thatformulation (1) and the solution algorithm proposed are general and applicable beyond the context of theFTCP, which we use solely as an example. Decision problems under endogenous uncertainty may, in fact,arise in several business context, some of which are mentioned in Section 1. To provide another practicalexample, consider an agribusiness making periodic production planning decisions for a number of differentcrops while demand is uncertain. Product substitution is a common phenomenon in the production of crops.In fact, a customer, say a farmer, may view multiple crops as suitable for their farming needs, see, e.g.,Bansal and Dyer (2020). Thus, the uncertain demand for a crop may depend in part on the portfolio ofcrops an agribusiness offers for sale, yielding a decision problem under endogenous uncertainty.In Section 4.1 we describe the FTCP in more details and formulate it as an MSPEU. In Section 4.2 weillustrate an efficient procedure to obtain big- M values. In Section 4.3 we introduce the problem instancesand finally in Section 4.4 we present and discuss the results. The data of the problem instances is avail-able online at https://github.com/GioPan/instancesFTCPwithEndogenousUncertainty for the benefitof future research. The FTCP is the problem of composing a football team by purchasing and selling professional footballplayers. A complete description of the problem can be found in Pantuso (2017). In what follows, we report8he basic elements necessary for this case study. The decision problem can be described as follows. At everystage, i.e., transfer market window (TMW), professional football clubs can renew their teams by purchasingplayers from other clubs or selling available players to other clubs. In order to participate in national andinternational competitions, professional clubs must compose a team made of a fixed number of players. Inaddition, the coach of the team requires players with different roles (e.g., defenders or mid-fielders) and skills.Club managers are often given a budget to spend in the transfer market and are typically allowedto reinvest the revenue from the sale of players. The current market value of football players is known,while the future value is stochastic as it depends on a number of random events such as injuries, fitness,motivation and ultimately luck. Therefore, at every TMW, football clubs make decisions in conditions ofuncertainty with the scope of maximizing the expected value of the team. Furthermore, the market valuesof the players in the same team are strongly correlated. This generates a multistage stochastic programwith endogenous uncertainty, since the decision of hiring or selling a player will change the correlations ofthe joint value distribution of the players considered. That is, hiring a football player will make their valuestrongly correlated with the value of the other players in the team, while selling or not hiring a player willmake their value uncorrelated with the value of the players in the team.In order to model the FTCP we assume the club is evaluating a finite number of alternative teamcompositions . Every team composition consists of the required number of players and contains the necessarymix of skills. Let I be the set of possible team compositions, N the set of nodes in the MDST describingthe underlying uncertainty, N in the set of child nodes of node n ∈ N if team composition i ∈ I is chosenat node n . Note in fact that the selection of a team composition determines the correlations between thevalues of the players in the instance, and thus the probability distribution. Let N L be the set of leaf nodes,i.e., the nodes at the last decision stage, and t ( n ) ∈ { , . . . , T } the stage corresponding to node n .Let V in be the value of team composition i ∈ I at node n (i.e., the sum of the values of the playersin team composition i ), let C Sin be the cost of the salary of the players in team composition i at node n , C Tjn the cost of transitioning from team composition i to team composition j at node n , i.e., the cost of thetransfer fees for the players bought minus the revenue for the transfer fees of the players sold. Furthermore,let B be the budget available to the club for the transfer market at node n and E n a stochastic extra-budgetconditional on events such as sport successes. Let X i be equal to one if team composition i ∈ I is the initialteam composition, M in a suitably large constant for each i and n (see Section 4.2), and ρ a discount rate.Let δ in be a binary decision variable which is equal 1 if team composition i ∈ I is chosen at node n ∈ N , x ijn a binary decision variable which is equal 1 if the club transitions from team composition i ∈ I to teamcomposition j ∈ J at node n ∈ N and, φ in the expected net value of the team at the children of node n ifcomposition i is chosen and, finally, θ n the expected net value of the team at the child nodes of node n ∈ N .The FTCP with endogenous uncertainty is hence:max (cid:88) i ∈I ( V i δ i − C Si δ i − (cid:88) j ∈I C Tij x ij ) + θ (6a)s.t. (cid:88) i ∈I δ in = 1 n ∈ N , (6b) δ i − X i + (cid:88) j ∈I : j (cid:54) = i ( x ij − x ji ) = 0 i ∈ I , (6c) δ in − δ i,a ( n ) + (cid:88) j ∈I : j (cid:54) = i ( x ijn − x jin ) = 0 n ∈ N \ { } , i ∈ I , (6d) (cid:88) i ∈I (cid:88) j (cid:54) = i ∈I x ijn ≤ n ∈ N , (6e)9 i ∈I (cid:88) j ∈I C Tijn x ijn ≤ B + E n n ∈ N , (6f) θ n ≤ φ in + M in (1 − δ in ) i ∈ I , n ∈ N , (6g) φ in = (cid:88) m ∈N in
11 + ρ t ( m ) π m θ m + (cid:88) j ∈I (cid:32) V jm δ jm − C Sjm δ jm − (cid:88) k ∈I C Tjkm x jkm (cid:33) i ∈ I , n ∈ N \ N L , (6h) θ n = 11 + ρ t ( n )+1 (cid:88) i ∈I V in δ in n ∈ N L , (6i) δ in ∈ { , } i ∈ I , n ∈ N , (6j) x ijn ∈ { , } i, j ∈ I , n ∈ N , (6k) θ n ∈ R n ∈ N . (6l)Objective function (6a) represents the sum of the value of the team composition chosen here-and-now,minus the expenses for salaries and, if any, the transition cost from the initial team composition. In additionit takes into account the expected value of the team at future nodes. Constraints (6b) ensure that onlyone team composition is chosen at each decision node, while constraints (6c)-(6d) state that either the clubholds the same team composition as in the previous season or a new team composition is chosen at node 0and at the rest of the nodes, respectively. Constraints (6e) state that at most one team change can be doneat each node. Constraints (6f) ensure that the net expenses for transitioning from a team composition toanother (i.e., the money spent in the transfer market) do not exceed the available budget. Constraints (6g)state that the future net expected value of the team depends on the team composition chosen. Constraints(6h) sets the expected value at the children of node n if composition i is chosen, while constraints (6i) setthe final value of θ n to the value of the team chosen at the leaf nodes. This corresponds to a sunset valuewhich accounts the termination of an infinite horizon problem. Notice that the final value of the team isdiscounted as a future value. Finally, constraints (6j)-(6l) define the domain for the decision variables. We illustrate a fast method to set the M in values in (6). The method proposed resulted, on this special case,significantly faster than the general method described in Appendix C.We start by observing that, for all i ∈ I and n ∈ N we are looking for a value M in such that: θ n − φ in ≤ M in and that, based on constraints (6g), we have θ n ≤ max i ∈I φ in consequently, we need to set M in ≥ max i ∈I φ in − φ in thus, under the very mild assumption that φ in ≥ i ∈ I a valid M in is M in = M n = max i ∈I φ in By expanding φ in according to constraints (6h) and (6i) we can set10 M in = max i ∈I ρ t ( n )+1 V in for all i ∈ I if n ∈ N L • M in = max i ∈I (cid:40) (cid:80) m ∈N in ρ t ( m ) π m (cid:18) M m + max j ∈I (cid:8) V jm (cid:9) (cid:19)(cid:41) for all i ∈ I if n ∈ N \ N L . Instances are generated from the case studies on the FTCP presented by Pantuso (2017) and are available on-line at https://github.com/GioPan/instancesFTCPwithEndogenousUncertainty . The instances consistof the 20 teams competing in the English Premier League 2013/2014 and dealing with the transfer marketof summer 2014. Each team is characterized by a budget, a list of players (made of the current and targetplayers), and a number of randomly generated team compositions I complying with the regulations of theleague and with the coach’s specifics. In turn, each player is characterized by age, role, current market value,salary, selling and purchase price. The value and cost of salaries of a team composition are calculated as thesum of values and salaries, respectively, of the players contained. Similarly, transition costs are calculatedas the revenue for the players sold, minus the cost of the players bought.Future player values are stochastic and modeled by means of the regression equation available in Pantuso(2017). The joint probability distribution of the market values for all the players forms a multivariatenormal distribution. In addition, we set a 0 . E n is calculated as a percentage of the deterministicbudget B . The percentage is the same as percentage increase of team value from the parent node, thatis max (cid:8) , ∗ ( V in − V i,a ( n ) ) /V i,a ( n ) (cid:9) . That is, an increase in the team value yields an increase in thespending capability of the team.The instances represent a four-stage horizon and are identified by (i) the focal team, where Team ∈ { ARS, ASV, CAR, CHE, CRP, EVE, FUL, HUL, LIV, MAC, MAU, NEC, NOR, SOU, STO, SUN, SWA,TOT, WBA, WHU } , (ii) the number of team composition |I| ∈ { , , } , and (iii) the number of samplesin each distribution S ∈ { , } , yielding 120 numerically different instances with up to approximately 400thousand binary variables. We implemented our algorithm in Python 2.7 using Cplex 12.8 for solving the subproblems. Cplex 12.8has also been used to solve the full problems. All experiments have been run on a server with 64 doubleAMD Opteron 6380 processors and 251 GB RAM. Tables 1 to 3 report the results for different cardinalities |I| and different number of samples S for all the focal teams. The computation times do not include thetime required to find big- M values as illustrated in Section 4.2. An account of these times is reportedin Appendix D. Note that big- M values are required both when solving the full problem and when usingour algorithm and in both cases are calculated beforehand. Thus, their computation time does not have animpact on the comparison between the two solution methods. Since the algorithm entails iteratively buildingand solving mixed-integer programs, the times reported include the time required for building the models.Table 1 reports the results with |I| = 3, generating instances with up to approximately 32000 binaryvariables. For these instances it can be noticed that Cplex performs much better than our method on thesmaller test cases (with S = 4 samples) while our algorithm is more competitive on the instances with S = 5.Altogether, Cplex solves the problems approximately 25% faster than our method.11able 1: Results with |I| = 3. S indicates the number of realizations describing each distribution. τ (Cpx) indicate the elapsed time when using the algorithm (Cplex). Obj. (Cpx)indicates the objective value obtained using the algorithm (Cplex). ∆ τ is calculated as 100( τ - τ Cpx)/ τ Cpx. Team S τ Cpx [sec] Obj. Cpx τ [sec] Obj. ∆ τ [%]ARS 4 24505 16965 22620 86.924 1487.985 178.413 1487.985 105.252ASV 4 24505 16965 22620 107.462 858.189 228.801 858.189 112.914CAR 4 24505 16965 22620 112.806 642.253 282.521 642.253 150.448CHE 4 24505 16965 22620 144.798 2604.929 297.918 2604.929 105.748CRP 4 24505 16965 22620 105.784 471.838 232.649 471.838 119.929EVE 4 24505 16965 22620 198.908 1216.299 189.362 1216.299 -4.799FUL 4 24505 16965 22620 116.068 694.942 187.200 694.942 61.284HUL 4 24505 16965 22620 174.406 571.997 215.792 571.997 23.729LIV 4 24505 16965 22620 457.376 1371.466 196.851 1371.466 -56.961MAC 4 24505 16965 22620 165.375 1734.123 166.470 1734.123 0.662MAU 4 24505 16965 22620 2036.260 2474.198 231.169 2474.198 -88.647NEC 4 24505 16965 22620 95.849 974.510 213.350 974.510 122.589NOR 4 24505 16965 22620 69.662 486.204 161.434 486.204 131.739SOU 4 24505 16965 22620 153.804 640.762 143.529 640.762 -6.681STO 4 24505 16965 22620 79.237 674.964 165.104 674.964 108.366SUN 4 24505 16965 22620 229.976 837.582 202.287 837.582 -12.040SWA 4 24505 16965 22620 90.615 659.840 228.751 659.840 152.442TOT 4 24505 16965 22620 193.830 1383.202 171.743 1383.202 -11.395WBA 4 24505 16965 22620 223.691 570.527 179.085 570.527 -19.941WHU 4 24505 16965 22620 259.587 785.027 202.723 785.027 -21.905ARS 5 47008 32544 43392 189.159 1508.956 403.177 1508.956 113.142ASV 5 47008 32544 43392 486.235 801.662 400.098 801.662 -17.715CAR 5 47008 32544 43392 223.908 661.053 389.787 661.053 74.083CHE 5 47008 32544 43392 593.065 2667.682 628.685 2667.682 6.006CRP 5 47008 32544 43392 222.443 478.977 483.007 478.977 117.138EVE 5 47008 32544 43392 726.302 1134.583 387.368 1134.583 -46.666FUL 5 47008 32544 43392 275.205 717.025 377.999 717.025 37.352HUL 5 47008 32544 43392 2785.695 494.981 446.281 494.981 -83.980LIV 5 47008 32544 43392 7386.953 1403.374 391.775 1403.374 -94.696MAC 5 47008 32544 43392 349.382 1769.971 306.545 1769.971 -12.261MAU 5 47008 32544 43392 7755.757 2260.280 488.009 2260.280 -93.708NEC 5 47008 32544 43392 481.964 1014.031 440.037 1014.031 -8.699NOR 5 47008 32544 43392 356.411 499.601 317.165 499.601 -11.012SOU 5 47008 32544 43392 148.499 656.613 322.197 656.613 116.968STO 5 47008 32544 43392 391.803 689.878 359.498 689.878 -8.245SUN 5 47008 32544 43392 481.092 854.091 406.024 854.091 -15.604SWA 5 47008 32544 43392 444.930 681.336 386.791 681.336 -13.067TOT 5 47008 32544 43392 418.714 1409.204 400.087 1409.204 -4.449WBA 5 47008 32544 43392 299.717 522.849 367.367 522.849 22.571WHU 5 47008 32544 43392 654.066 678.902 426.684 678.902 -34.764Avg 744.342 305.093 25.378By increasing the number of team compositions we generate larger MDSTs, and the correspondingMSPEUs also become larger. Table 2 reports the results with |I| = 4 generating a MDST with four possibledistributions at each decision node and problems with up to approximately 135 thousand binary variables.In this case, our algorithm performs significantly better that Cplex on almost all instances, solving the12roblems 32.8% faster.Table 2: Results with |I| = 4. S indicates the number of realizations describing each distribution. τ (Cpx) indicate the elapsed time when using the algorithm (Cplex). Obj. (Cpx)indicates the objective value obtained using the algorithm (Cplex). ∆ τ is calculated as 100( τ - τ Cpx)/ τ Cpx.Case S τ Cpx [sec] Obj. Cpx τ [sec] Obj. ∆ τ [%]ARS 4 91749 69904 65535 469.440 1409.074 585.920 1409.074 24.812ASV 4 91749 69904 65535 560.788 916.248 453.329 916.248 -19.162CAR 4 91749 69904 65535 551.781 574.859 524.636 574.859 -4.920CHE 4 91749 69904 65535 756.516 2476.158 729.832 2476.158 -3.527CRP 4 91749 69904 65535 558.236 431.065 632.984 431.065 13.390EVE 4 91749 69904 65535 7656.373 1282.401 535.382 1282.401 -93.007FUL 4 91749 69904 65535 932.283 711.560 488.729 711.560 -47.577HUL 4 91749 69904 65535 7746.679 559.913 608.927 559.913 -92.140LIV 4 91749 69904 65535 7627.856 1452.961 545.440 1452.961 -92.849MAC 4 91749 69904 65535 888.410 1657.326 371.674 1657.326 -58.164MAU 4 91749 69904 65535 7672.755 2489.203 701.230 2489.203 -90.861NEC 4 91749 69904 65535 545.592 908.271 518.564 908.271 -4.954NOR 4 91749 69904 65535 393.393 443.482 371.138 443.482 -5.657SOU 4 91749 69904 65535 368.125 603.976 342.766 603.976 -6.889STO 4 91749 69904 65535 453.101 633.148 396.288 633.148 -12.539SUN 4 91749 69904 65535 529.445 785.688 513.954 785.688 -2.926SWA 4 91749 69904 65535 469.695 607.176 440.718 607.176 -6.169TOT 4 91749 69904 65535 639.796 1295.453 421.155 1295.453 -34.174WBA 4 91749 69904 65535 3945.793 633.948 473.058 633.948 -88.011WHU 4 91749 69904 65535 8498.951 828.664 601.910 828.664 -92.918ARS 5 176841 134736 126315 1210.558 1440.502 1340.235 1440.502 10.712ASV 5 176841 134736 126315 1946.890 760.043 1381.227 760.043 -29.055CAR 5 176841 134736 126315 1357.787 597.870 1329.530 597.870 -2.081CHE 5 176841 134736 126315 1706.027 2602.808 1896.954 2602.808 11.191CRP 5 176841 134736 126315 1365.008 467.102 1489.021 467.102 9.085EVE 5 176841 134736 126315 8505.074 1232.924 1260.209 1232.924 -85.183FUL 5 176841 134736 126315 7808.254 699.597 1388.999 699.597 -82.211HUL 5 176841 134736 126315 8642.826 553.807 1300.262 553.807 -84.956LIV 5 176841 134736 126315 8326.633 1460.389 1234.644 1460.389 -85.172MAC 5 176841 134736 126315 1020.989 1678.841 1149.133 1678.841 12.551MAU 5 176841 134736 126315 8397.024 2285.317 1405.108 2285.317 -83.267NEC 5 176841 134736 126315 1623.349 1062.903 1427.509 1062.903 -12.064NOR 5 176841 134736 126315 1041.199 480.087 1085.893 480.087 4.293SOU 5 176841 134736 126315 1000.712 616.957 1020.687 616.957 1.996STO 5 176841 134736 126315 1111.253 663.486 1059.614 663.486 -4.647SUN 5 176841 134736 126315 1480.353 891.809 1609.855 891.809 8.748SWA 5 176841 134736 126315 1197.202 623.004 1596.017 623.004 33.312TOT 5 176841 134736 126315 2560.338 1441.642 1264.625 1441.642 -50.607WBA 5 176841 134736 126315 8304.973 601.984 1152.420 601.984 -86.124WHU 5 176841 134736 126315 8432.614 755.708 1416.893 755.708 -83.197Avg 3207.601 926.661 -32.873This pattern illustrates that our algorithm scales significantly better than Cplex, and is confirmed also byour results on the instances with |I| = 5, the largest we tested. These instances generate problems with up13o approximately 400 thousand binary variables and the corresponding results are provided in Table 3. It ispossible to notice that our algorithm outperforms Cplex on almost all instances. Altogether, our algorithmis able to solve the same instances 61.1% faster than Cplex. Only on two instances (CAR and CHE, with S = 4) Cplex performs better than our algorithm. Similar cases are to be expected since our algorithmentails solving several subproblems which are mixed-integer programs. Therefore, it is possible that somenumerically difficult subproblems are found that slow down the entire algorithm (in our case we do notsolve problems in parallel). However, despite that, the average performance of our algorithm on the largestinstance is by far better.These results illustrate that MSPEU easily generate very large optimization problems. The size of theproblems, in our instances, increases by approximately three to five times by increasing the number ofpossible distributions. However, the algorithm we provide scales better than the solver Cplex. Particularly,when using Cplex, the average solution time increases by approximately 330% when increasing |I| from 3to 4, and by approximately 175% when increasing |I| from 4 to 5. With our algorithm the average solutiontime increases by approximately 203% when increasing |I| from 3 to 4, and by approximately 158% whenincreasing |I| from 4 to 5.Table 3: Results with |I| = 5. S indicates the number of realizations describing each distribution. τ (Cpx) indicate the elapsed time when using the algorithm (Cplex). Obj. (Cpx)indicates the objective value obtained using the algorithm (Cplex). ∆ τ is calculated as 100( τ - τ Cpx)/ τ Cpx.Case S τ Cpx [sec] Obj. Cpx τ [sec] Obj. ∆ τ [%]ARS 4 261051 210525 151578 4148.394 1412.322 2487.630 1412.322 -40.034ASV 4 261051 210525 151578 9128.747 888.394 2360.369 888.394 -74.144CAR 4 261051 210525 151578 2324.632 559.909 2951.872 559.909 26.982CHE 4 261051 210525 151578 2941.786 2351.709 4074.808 2351.709 38.515CRP 4 261051 210525 151578 4873.336 403.008 2757.925 403.008 -43.408EVE 4 261051 210525 151578 9240.398 1318.359 2367.225 1328.380 -74.382FUL 4 261051 210525 151578 9385.279 676.308 2775.651 676.308 -70.425HUL 4 261051 210525 151578 9616.179 496.264 2430.993 600.252 -74.720LIV 4 261051 210525 151578 9892.436 1584.601 1522.159 1584.601 -84.613MAC 4 261051 210525 151578 1748.831 1589.133 1070.475 1589.133 -38.789MAU 4 261051 210525 151578 9294.710 2529.355 1935.253 2533.948 -79.179NEC 4 261051 210525 151578 6984.354 1202.156 1502.823 1202.156 -78.483NOR 4 261051 210525 151578 3516.805 420.312 1038.191 420.312 -70.479SOU 4 261051 210525 151578 3593.291 599.703 947.766 599.703 -73.624STO 4 261051 210525 151578 1966.818 632.875 1182.410 632.875 -39.882SUN 4 261051 210525 151578 4767.951 731.767 1460.145 731.767 -69.376SWA 4 261051 210525 151578 4397.444 576.314 1280.492 576.314 -70.881TOT 4 261051 210525 151578 3894.520 1392.754 1169.367 1392.754 -69.974WBA 4 261051 210525 151578 11192.050 590.701 1298.701 661.566 -88.396WHU 4 261051 210525 151578 9314.423 600.960 1755.932 768.724 -81.148ARS 5 504556 406900 292968 5985.321 1453.229 3231.294 1453.229 -46.013ASV 5 504556 406900 292968 13100.017 969.935 2557.610 969.935 -80.476CAR 5 504556 406900 292968 12990.803 612.171 3115.356 612.171 -76.019CHE 5 504556 406900 292968 15815.020 2608.034 4221.096 2608.034 -73.310CRP 5 504556 406900 292968 13160.004 476.955 3668.813 476.955 -72.121EVE 5 504556 406900 292968 18368.885 1289.379 2913.376 1452.827 -84.140FUL 5 504556 406900 292968 13385.821 746.554 2811.966 755.759 -78.993HUL 5 504556 406900 292968 12995.828 516.157 3335.790 655.447 -74.332LIV 5 504556 406900 292968 12808.767 1473.064 2854.447 1706.759 -77.715MAC 5 504556 406900 292968 5341.258 1681.019 2050.793 1681.019 -61.605MAU 5 504556 406900 292968 13282.279 2690.920 3722.067 2712.864 -71.97714EC 5 504556 406900 292968 13368.708 1365.929 2971.268 1365.929 -77.774NOR 5 504556 406900 292968 5352.415 462.292 2111.632 462.292 -60.548SOU 5 504556 406900 292968 5347.118 624.205 1897.896 624.205 -64.506STO 5 504556 406900 292968 5679.705 660.999 2308.348 660.999 -59.358SUN 5 504556 406900 292968 13483.620 927.387 2961.979 927.387 -78.033SWA 5 504556 406900 292968 11568.065 631.495 2500.653 631.495 -78.383TOT 5 504556 406900 292968 12955.368 1490.643 2323.134 1490.643 -82.068WBA 5 504556 406900 292968 12852.950 709.452 2696.335 718.245 -79.022WHU 5 504556 406900 292968 13408.572 843.851 3371.308 844.027 -74.857Avg 8836.822 2399.883 -65.192 This paper introduced 1) a novel scenario tree structure and 2) a node formulation for multistage stochasticprograms with endogenous uncertainty, as well as 3) a solution algorithm for a special case. A computa-tional study shows that while, as expected, problems with endogenous uncertainty tend to generate largeoptimization problems, all our instances where solvable by Cplex to optimality in at most approximately 5hours. Furthermore, our algorithm outperformed Cplex on the medium and large instances and showed thatit scales well with the size of the problem.Despite the encouraging results obtained in our study, solving multistage stochastic programs with en-dogenous uncertainty remains, in general, a challenging task. Our algorithm requires an explicit scenario treestructure, and solves a number of problems which grows linearly with the number of nodes. However, thenumber of nodes in a scenario tree grows exponentially with the number of stages and the treatment of caseswith more than a handful of stages may soon become prohibitive. New approaches in the spirit of Pereiraand Pinto (1991); Zou et al. (2019), based on progressive approximations of future stages, may be provenmore scalable. Furthermore, our models employ so called big- M constants. Poorly chosen big- M values,e.g., by trial-and-error, may become problematic. It is well known that they may create numerical difficultieswhen solving mixed-integer programs. In addition, as discussed in Pineda and Morales (2019), they maylead to highly sub-optimal solutions. The authors use a simple bilevel programming problem (which can bereformulated as a mixed-integer program that uses big- M s) to show how a poorly designed trial-and-errorprocedure may generate the false belief that the solution to the reformulation is indeed optimal for the orig-inal bilevel problem. Based on these evidences we also advocate caution and the use of more sophisticatedprocedures for setting big- M values. The procedure in Appendix C goes in this direction. Finally, casesmore general than the special one treated Section 3 remain to be addressed. References
S. Ahmed.
Strategic planning under uncertainty: Stochastic integer programming approaches . Phd, Universityof Illinois at Urbana-Champaign, 2000.R. M. Apap and I. E. Grossmann. Models and computational strategies for multistage stochastic program-ming under endogenous and exogenous uncertainties.
Computers & Chemical Engineering , 103:233–274,2016. ISSN 0098-1354.S. Bansal and J. S. Dyer. Planning for end-user substitution in agribusiness.
Operations Research , 68(4):1000–1019, 2020.J. R. Birge and F. Louveaux.
Introduction to stochastic programming . Springer, New York, 1997.M. Colvin and C. T. Maravelias. A stochastic programming approach for clinical trial planning in new drugdevelopment.
Computers & Chemical Engineering , 32(11):2626–2642, 2008.M. Colvin and C. T. Maravelias. Modeling methods and a branch and cut algorithm for pharmaceuticalclinical trial planning using stochastic programming.
European Journal of Operational Research , 203(1):205–215, 2010. 15. da Costa Flach.
Stochastic Programming with Endogenous Uncertainty: An Application in HumanitarianLogistics . PhD thesis, PUC-Rio, 2010.L. F. Escudero, M. A. Gar´ın, J. F. Monge, and A. Unzueta. On preparedness resource allocation planningfor natural disaster relief under endogenous uncertainty with time-consistent risk-averse management.
Computers & Operations Research , 98:84–102, 2018.V. Goel and I. E. Grossmann. A stochastic programming approach to planning of offshore gas field develop-ments under uncertainty in reserves.
Computers & Chemical Engineering , 28:1409–1429, 2004.V. Goel and I. E. Grossmann. A Class of stochastic programs with decision dependent uncertainty.
Mathe-matical Programming , 108:355–394, 2006.V. Gupta and I. E. Grossmann. Solution strategies for multistage stochastic programming with endogenousuncertainties.
Computers & Chemical Engineering , 35:2235–2247, 2011.H. Held and D. L. Woodruff. Heuristics for Multi-Stage Interdiction of Stochastic Networks.
Journal ofHeuristics , 11(5-6):483–500, 2005.L. Hellemo.
Managing Uncertainty in Design and Operation of Natural Gas Infrastructure . Phd, NorwegianUniversity of Science and Technology, 2016.L. Hellemo, P. I. Barton, and A. Tomasgard. Decision-dependent probabilities in stochastic programs withrecourse.
Computational Management Science , 15(3):369–395, 2018.F. Hooshmand and S. A. MirHassani. Efficient constraint reduction in multistage stochastic programmingproblems with endogenous uncertainty.
Optimization Methods and Software , 31:359–376, 2016.F. Hooshmand and S. A. MirHassani. Reduction of nonanticipativity constraints in multistage stochasticprogramming problems with endogenous and exogenous uncertainty.
Mathematical Methods of OperationsResearch , 87(1):1–18, 2018.T. W. Jonsbr˚aten, R. J.-B. Wets, and D. L. Woodruff. A class of stochastic programs withdecision dependentrandom elements.
Annals of Operations Research , 82:83–106, 1998.P. Kall and S. W. Wallace.
Stochastic programming . John Wiley and Sons Ltd, Chichester, 1994.M. Kaut, K. T. Midthun, A. S. Werner, A. Tomasgard, L. Hellemo, and M. Fodstad. Multi-horizon stochasticprogramming.
Computational Management Science , 11(1-2):179–193, 2014.M. Laumanns, S. Prestwich, and B. Kawas. Distribution shaping and scenario bundling for stochasticprograms with endogenous uncertainty. In J. L. Higle, W. R¨omisch, and S. Sen, editors,
StochasticProgramming E-Print Series , Stochastic Programming E-Print Series. Institut f¨ur Mathematik, 2014.L. Mercier and P. Van Hentenryck. An anytime multistep anticipatory algorithm for online stochasticcombinatorial optimization.
Annals of Operations Research , 184(1):233–271, 2011.G. Pantuso. The Football Team Composition Problem: A Stochastic Programming approach.
Journal ofQuantitative Analysis in Sports , 13(3):113–129, 2017.G. Pantuso and L. M. Hvattum. Maximizing performance with an eye on the finances: a chance-constrainedmodel for football transfer market decisions.
TOP , pages 1–29, 2020.S. Peeta, F. S. Salman, D. Gunnec, and K. Viswanath. Pre-disaster investment decisions for strengtheninga highway network.
Computers & Operations Research , 37:1708–1719, 2010.M. V. Pereira and L. M. Pinto. Multi-stage stochastic optimization applied to energy planning.
Mathematicalprogramming , 52(1-3):359–375, 1991.S. Pineda and J. M. Morales. Solving linear bilevel problems using big-ms: Not all that glitters is gold.
IEEE Transactions on Power Systems , 34(3):2469–2471, 2019. doi: 10.1109/TPWRS.2019.2892607.16. Tarhan and I. E. Grossmann. A multistage stochastic programming approach with strategies for un-certainty reduction in the synthesis of process networks with uncertain yields.
Computers & ChemicalEngineering , 32:766–788, 2008.B. Tarhan, I. E. Grossmann, and V. Goel. Stochastic Programming Approach for the Planning of Offshore Oilor Gas Field Infrastructure under Decision-Dependent Uncertainty.
Industrial & Engineering ChemistryResearch , 48:3078–3097, 2009.B. Tarhan, I. E. Grossmann, and V. Goel. Computational strategies for non-convex multistage MINLPmodels with decision-dependent uncertainty and gradual uncertainty resolution.
Annals of OperationsResearch , 203(1):141–166, 2013.K. Tong, Y. Feng, and G. Rong. Planning under Demand and Yield Uncertainties in an Oil Supply Chain.
Industrial & Engineering Chemistry Research , 51(2):814–834, 2012.K. Viswanath, P. Srinivas, and S. F. Sibel. Investing in the Links of a Stochastic Network to Minimize Ex-pected Shortest Path Length. 2004. URL .J. Zou, S. Ahmed, and X. A. Sun. Stochastic dual dynamic integer programming.
Mathematical Programming ,175(1-2):461–502, 2019.
A Notation table
Table 4: Notation of problem (1).Sets { , . . . , T } Set of decision stages N Set of nodes in the multi-distribution scenario tree N t ⊆ N Set of nodes at stage t D n Set of possible distributions applicable at node n N nd ⊆ N Set of child nodes of node n if distribution d ∈ D n is enforced X t Domain of the decision variables at stage t Parameters t ( n ) Stage of node na ( n ) Parent node of node nπ n Probability of node nr n ∈ R N t ( n ) Coefficients of decision variables x n in the objective function q nd ∈ R Coefficient of decision variable δ nd in the objective function A n ∈ R M t ( n ) × N t ( n ) Coefficients of variables x n in the constraints that connect x n and δ nd decisions B nd ∈ R M t ( n ) × Coefficients of variable δ nd in the constraints that connect x n and δ nd decisions C a ( n ) ∈ R M t ( n ) × N t ( a ( n )) Coefficients of variables x a ( n ) in the constraints that connect x n and δ nd decisions D a ( n ) ,d ∈ R M t ( n ) × Coefficients of variable δ a ( n ) ,d in the constraints that connect x n and δ nd decisions h n ∈ R M t ( n ) Right-hand-side coefficients of the constraints that connect x n and δ nd decisionsΘ n ∈ R Terminal value of the decisions following node n ∈ N T Variables x n ∈ R N t ( n ) Decisions made at node nδ nd ∈ { , } Decision on whether to apply probability distribution d at node nθ n ∈ R Expected value of the decisions made at the nodes descending from n A big- M reformulation In this appendix a big- M reformulation that linearizes model (1) is introduced. In addition to the notationintroduced in Section 2, let M nd ∈ R be a suitably high constant. The linearized EUMSP is thusmax r T x + (cid:88) d ∈D q d δ d + θ (7a)s.t. (cid:88) d ∈D n δ nd = 1 n ∈ N , (7b) A n x n + (cid:88) d ∈D n B nd δ nd + C a ( n ) x a ( n ) + (cid:88) d ∈D a ( n ) D a ( n ) ,d δ a ( n ) ,d = h n n ∈ N , (7c) θ n ≤ (cid:88) m ∈N nd π m ( r Tm x m + (cid:88) d ∈D m q md δ md + θ m ) + M nd (1 − δ nd ) n ∈ N \ N T , d ∈ D n (7d) θ n = Θ n n ∈ N T , (7e) x n ∈ X t ( n ) n ∈ N , (7f) δ nd ∈ { , } n ∈ N , d ∈ D n , (7g) θ n ∈ R n ∈ N . (7h)Note, particularly, that constraints (7d) are equivalent to (1d). Consider a given node n , other than a leafnode. Observe that only for one distribution d there will be a δ nd which takes value one at n (see (7b)). Forthe same n and for the same d , the second term on the right-hand-side of (7d) will be zero (i.e., the big- M will not be enforced), and the resulting right-hand-side will be the most binding among the |D n | constraintsfor node n . Since we are maximizing, at optimality θ n will take value of the expectation according to thedistribution d for which δ nd = 1, as it happens in model (1). C Finding big- M values An efficient implementation of model (7) requires tight big- M values. Observe that, for n ∈ N \ N T and d ∈ D n , constant M nd must be a valid upper bound for constraints (7d), that is: θ n − (cid:88) m ∈N nd π m ( r Tm x m + (cid:88) d ∈D m q md δ md + θ m ) ≤ M nd Let us introduce φ nd to represent the expectation at the children of node n for distribution d , that is: φ nd = (cid:88) m ∈N nd π m ( r Tm x m + (cid:88) d ∈D m q md δ md + θ m )Consider the numerical example shown in Table 5 for a given node n and three possible distributions d . Thetable reports the highest and lowest values the expectation at the following stage can take for each possibledistribution, and the corresponding value of θ n should a specific distribution be chosen. When choosing M n,d , notice that the maximum value θ n can reach for other distributions is 9, and that the least value φ n,d can reach is 5. Therefore, we can set M n,d = 9 −
5. In fact, if distribution d is chosen, θ n will beat most 9 and φ n,d at least 5, thus adding 4 to φ n,d will ensure that θ n is correctly set to 9. Similarly, wechoose M n,d = 6 as the highest value θ n can take if d is not selected is 10, while the least value of φ n,d is4. Finally, with a similar reasoning we can set M n,d = 7.From the example in Table 5 we understand that finding values for M nd amounts to finding highest valuesfor φ nd and differences θ n − φ nd . In what follows we illustrate how these values can be found for t = T − t = T − , . . . , C.1 Big- M values for stages t = T − We start at the second-last stage, t = T −
1. Our task is that of finding, for each node ¯ n ∈ N T − and foreach distribution ¯ d ∈ D ¯ n , a constant M ¯ n ¯ d which is slightly higher than the highest difference θ ¯ n − φ ¯ n ¯ d , where18able 5: Numerical example for the calculation of constants M nd φ nd d θ n Max Min M nd d
10 10 5 9 − d − d − φ ¯ n ¯ d = (cid:88) m ∈N ¯ n ¯ d π m ( r Tm x m + (cid:88) d ∈D m q md δ md + θ m )Now, the highest difference can be found solving the following optimization problem: M ∗ ¯ n ¯ d = max θ ¯ n − (cid:88) m ∈N ¯ n ¯ d π m ( r Tm x m + (cid:88) k ∈D m q mk δ mk + Θ m ) (8a)s.t. (cid:88) d ∈D n δ nd = 1 n ∈ N , (8b) A n x n + (cid:88) d ∈D n B nd δ nd + C a ( n ) x a ( n ) + (cid:88) d ∈D a ( n ) D a ( n ) ,d δ a ( n ) ,d = h n n ∈ N , (8c) θ ¯ n ≤ Θ ∗ ¯ n ¯ d , (8d) x n ∈ X t ( n ) n ∈ N , (8e) δ nd ∈ { , } n ∈ N , d ∈ D n , (8f) δ ¯ n ¯ d = 0 (8g)Problem (8) consists of finding the feasible solution to problem (7) which yields the highest value for theleft-hand-side of constraint (7d) for ¯ n and ¯ d . The following two elements must be noted in (8). Constraints(7d) of the original problem, which determine the correct expectations at the stages before T −
1, are notincluded as they are irrelevant for stage T −
1. The second element to note is constraint (8d) which sets anupper bound θ ¯ n . This upper bound represents the highest value θ ¯ n can take for the distributions other than¯ d . This value can, in turn, be obtained solving optimization problems. The highest expectation for stage T ,given distribution d (cid:48) ∈ D ¯ n is the optimal value to problem (9):Φ ∗ ¯ nd (cid:48) = max (cid:88) m ∈N ¯ nd (cid:48) π m ( r Tm x m + (cid:88) k ∈D m q mk δ mk + Θ m ) (9a)s.t. (cid:88) d ∈D n δ nd = 1 n ∈ N , (9b) A n x n + (cid:88) d ∈D n B nd δ nd + C a ( n ) x a ( n ) + (cid:88) d ∈D a ( n ) D a ( n ) ,d δ a ( n ) ,d = h n n ∈ N , (9c) x n ∈ X t ( n ) n ∈ N , (9d) δ nd ∈ { , } n ∈ N , d ∈ D n , (9e) δ ¯ nd (cid:48) = 1 (9f)Therefore, when calculating M ¯ n ¯ d , the upper bound Θ ∗ ¯ n in (8d) is given by:Θ ∗ ¯ n ¯ d = max d ∈D ¯ n : d (cid:54) = ¯ d Φ ∗ ¯ nd Clearly, solving problems (8) and (9) amounts to solving integer programs of size comparable with theoriginal problem (7). However, the tightest Θ ∗ ¯ n ¯ d and M ¯ n ¯ d are not necessary, and higher values would stillprovide correct results. A suitable value for M ¯ n ¯ d can be obtained by solving any relaxation of problem (8),19ielding M R ¯ n ¯ d , and (9), yielding Φ R ¯ n ¯ d and in turn Θ R ¯ n ¯ d . As an example, one might solve the linear programmingrelaxation of problems (8) and (9) or, if the size of the problems is excessively high, one might choose torelax constraints (8c) and (9c) for some stages. Finally, since the procedure outlined might return negativevalues for some M nd , we set M nd = max { , M Rnd } to reduce, when possible, high big- M absolute values. Theprocedure is summarized in Algorithm 1. Algorithm 1
Algorithm for calculating M nd for T − Input: N , D n for n ∈ N , Θ n for n ∈ N T for Node ¯ n ∈ N T − do for Distribution d ∈ D ¯ n do Calculate Φ LP ¯ nd by solving the LP relaxation to problem (9) end for for Distribution ¯ d ∈ D ¯ n do In constraint (8d) set Θ ∗ ¯ n ¯ d = max d ∈D ¯ n : d (cid:54) = ¯ d Φ LP ¯ nd Calculate M R ¯ n ¯ d by solving a suitable relaxation to problem (8) Set M ¯ n ¯ d = max { , M R ¯ n ¯ d } end for end for return M nd for n ∈ N T − and d ∈ D n . C.2 Big- M values for stages t = T − , . . . , Once constants M nd are available for every n ∈ N T − and d ∈ D n we can proceed in a similar way tocalculate big- M s for stages T − , . . . ,
1. Given a stage ¯ t ∈ { T − , . . . , } , a node at that stage, ¯ n ∈ N ¯ t , anddistribution available at that node ¯ d ∈ D ¯ n , the tightest value of constant M ¯ n ¯ d , namely M ∗ ¯ n ¯ d , is M ∗ ¯ n ¯ d = max θ ¯ n − (cid:88) m ∈N ¯ n ¯ d π m ( r Tm x m + (cid:88) k ∈D m q mk δ mk + θ m ) (10a)s.t. (cid:88) d ∈D n δ nd = 1 n ∈ N , (10b) A n x n + (cid:88) d ∈D n B nd δ nd + C a ( n ) x a ( n ) + (cid:88) d ∈D a ( n ) D a ( n ) ,d δ a ( n ) ,d = h n n ∈ N , (10c) θ n ≥ (cid:88) m ∈N nd π m ( r Tm x m + (cid:88) k ∈D m q mk δ mk + θ m ) − M nd (1 − δ nd ) t = ¯ t + 1 , . . . , T − , n ∈ N t , d ∈ D n , (10d) θ n = Θ n n ∈ N T (10e) θ ¯ n ≤ Θ ∗ ¯ n ¯ d , (10f) x n ∈ X t ( n ) n ∈ N , (10g) δ nd ∈ { , } n ∈ N , d ∈ D n , (10h) δ ¯ n ¯ d = 0 (10i) Problem (10) consists of finding the feasible solution to problem (7) which yields the highest value of theleft-hand-side of constraint (7d) for node ¯ n and distribution ¯ d . Notice that, unlike in problem (8), problem(10) includes constraints (10d) which are necessary to ensure that θ n values are set to the lowest expectationfor all stages between ¯ t and T −
1. Note that constant M nd in constraints (10d) is also an upper bound to thequantity θ n − (cid:80) m ∈N nd π m ( r Tm x m + (cid:80) k ∈D m q mk δ mk + θ m ) and can be thus set to the quantities determinedat previous iterations. Also in this case, Θ ∗ ¯ n in (10f) represents the highest possible value θ ¯ n can take fordistributions other than ¯ d . The highest expectation for stage ¯ t + 1, given distribution d (cid:48) ∈ D ¯ n is the optimalvalue to problem (11) Φ ∗ ¯ nd (cid:48) = max (cid:88) m ∈N ¯ nd (cid:48) π m ( r Tm x m + (cid:88) k ∈D m q mk δ mk + θ m ) (11a)s.t. (cid:88) d ∈D n δ nd = 1 n ∈ N , (11b) n x n + (cid:88) d ∈D n B nd δ nd + C a ( n ) x a ( n ) + (cid:88) d ∈D a ( n ) D a ( n ) ,d δ a ( n ) ,d = h n n ∈ N , (11c) θ n ≤ (cid:88) m ∈N nd π m ( r Tm x m + (cid:88) k ∈D m q mk δ mk + θ m ) + M nd (1 − δ nd ) t = ¯ t + 1 , . . . , T − , n ∈ N t , d ∈ D n , (11d) θ n = Θ n n ∈ N T (11e) x n ∈ X t ( n ) n ∈ N , (11f) δ nd ∈ { , } n ∈ N , d ∈ D n , (11g) δ ¯ nd (cid:48) = 1 (11h) Therefore, the upper bound Θ ∗ ¯ n ¯ d in (10f) is given by:Θ ∗ ¯ n = max d ∈D ¯ n : d (cid:54) = ¯ d Φ ∗ ¯ nd Similarly to Section C.1, calculating the optimal M ∗ ¯ n ¯ d is cumbersome as well as not strictly necessary.Therefore, any computationally suitable relaxation to problems (10) and (11) can be adopted. The procedurefor obtaining constants M nd for stages T − , . . . , Algorithm 2
Algorithm for calculating M nd for ¯ t = T − , . . . , Input: N , D n for n ∈ N , Θ n for n ∈ N T , M nd for n ∈ N t , t = ¯ t + 1 , . . . , T − d ∈ D n . for Node ¯ n ∈ N ¯ t do for Distribution d ∈ D ¯ n do Calculate Φ LP ¯ nd by solving the LP relaxation to problem (11) end for for Distribution ¯ d ∈ D ¯ n do In constraint (10f) set Θ ∗ ¯ n ¯ d = max d ∈D ¯ n : d (cid:54) = ¯ d Φ LP ¯ nd Calculate M R ¯ n ¯ d by solving a relaxation to problem (10) Set M ¯ n ¯ d = max { , M R ¯ n ¯ d } end for end for return M nd for n ∈ N ¯ t and d ∈ D n . D Computation time for finding big- M values for the FTCP Table 6: Average elapsed time in seconds for the computations of big- M values using the procedure inSection 4.2. S indicates the number of realizations describing each distribution.Team S |I| = 3 |I| = 4 |I||I|