- Email: [email protected]
A Riccati-based primal interior point solver for multistage stochastic programming
European Journal of Operational Research 143 (2002) 452–461 www.elsevier.com/locate/dsw
A Riccati-based primal interior point solver for multistage stochastic programming €rgen Blomvall *, Per Olov Lindberg Jo Department of Mathematics, Link€opings Universitet, SE-58183 Link€oping, Sweden Received 11 January 2000; accepted 12 July 2001
Abstract We propose a new method for certain multistage stochastic programs with linear or nonlinear objective function, combining a primal interior point approach with a linear-quadratic control problem over the scenario tree. The latter problem, which is the direction finding problem for the barrier subproblem is solved through dynamic programming using Riccati equations. In this way we combine the low iteration count of interior point methods with an efficient solver for the subproblems. The computational results are promising. We have solved a financial problem with 1,000,000 scenarios, 15,777,740 variables and 16,888,850 constraints in 20 hours on a moderate computer. Ó 2002 Elsevier Science B.V. All rights reserved. Keywords: Dynamic programming; Finance; Interior point methods; Stochastic programming
1. Introduction Stochastic programming has many important applications (see e.g. Birge, 1997). It further easily leads to huge mathematical programs. It would thus be a natural candidate for the application of interior point methods. However, a direct application typically leads to dense systems of normal equations, essentially preventing this approach. Some authors have worked on circumventing these difficulties (Birge and Holmes, 1992) including work aiming at parallel computations (Vladimirou and Zenios, 1999). See also Birge (1997) and the references therein. *
Corresponding author. Tel.: +46-13-281406; fax: +46-13285770. E-mail address: [email protected] (J. Blomvall).
Others (e.g. Gondzio and Kouwenberg, 2001) have used interior point methods to solve the problems that appear in decomposition schemes for stochastic programming, such as L-shaped decompostion (Van Slyke and Wets, 1969). In the current paper we apply a primal interior point method to certain stochastic programs. Our main source of motivation comes from the financial area. Our work is related to, but independent of that of Steinbach (2000), whose approach is based on primal–dual interior point methods. Our approach is based on viewing the problem of computing a Newton step in the barrier subproblem as a linear-quadratic control problem (LQCP) over the scenario tree. For standard LQCPs there are standard techniques based on Riccati equations to compute closed form solutions (see e.g. Bertsekas, 1995). These techniques carry over
0377-2217/02/$ - see front matter Ó 2002 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 7 - 2 2 1 7 ( 0 2 ) 0 0 3 0 1 - 6
J. Blomvall, P.O. Lindberg / European Journal of Operational Research 143 (2002) 452–461
to our situation as we will show. Salinger and Rockafellar (1999) in applying operator splitting to stochastic programs also arrive at LQCPs over trees and observe that they can be treated by Riccati equations. In the rest of Section 1, we will give our model framework and state our model. Section 2 is devoted to deriving the algorithm. Our test problem, which motivated this study, is given in Section 3 together with some computational results. We are able to solve quite large problems (1,000,000 scenarios giving 15,777,740 nontrivial variables and 16,888,850 constraints) in reasonable time (20 hours) with rather limited computer resources. 1.1. Model framework Stochastic programs with discrete probability distribution can equivalently be stated as mathematical programs. These can be represented by time-staged trees, scenario trees (see Fig. 1), where each node represents a certain outcome of the random variables up to that stage. The stages (the depth from the root node) represent the evolution of time. Each node i 2 I has one predecessor
node, i , and a set, Iþ i , of immediate successor nodes in the next stage. (This is the same notation as in Salinger and Rockafellar (1999).) The leftmost node, f1g, in the tree is the root node. In the first time stage there is no uncertainty. The terminal nodes, T, represent the last time stage. The set Iþþ contains all the successors of node i, i.e. all i successors, up to and including the leaf nodes. 1.2. Model The notation above allows for convenient formulation of many stochastic programming problems, for example in finance which has been our main source of inspiration. Assume that we want to solve the following stochastic programming problem (P), partly in ‘‘control form’’: ðP Þ min xi ;ui
X
Ui ðxi ; ui Þ
Ei x i þ F i u i P e i ; ui 2
ð1aÞ
i2I
s:t: xi ¼ Ai xi þ Bi ui þ bi ; xi 2
Fig. 1. Notation for a multistage model.
453
½xli ; xui ; ½uli ; uui ;
i 2 I;
i 2 I n f1g;
ð1bÞ ð1cÞ
i 2 I;
ð1dÞ
i 2 I;
ð1eÞ
where x1 is a constant. We say that the problem is in control form, since the decision variables, are partitioned into state variables xi and control variables ui . For any values of the controls fui gi2I the states fxi gi2I are uniquely determined through state transitions equations (1b) (note that (1b) is of the same format as in a discrete time control problem). On top of (1b), we have local (node) constraints (1c) and simple bounds (1d)–(1e). The objective is a sum of convex functions of the node variables. The probability of each node is assumed to be included in Ui ð Þ. Many stochastic programs can be formulated as (P). In our main financial application, the ui are buying and selling decisions, whereas the xi are ‘‘inventories’’ of stocks, options, cash etc. Letting j j denote the number of elements of a set then we can define the following composite variables: xT ¼ ðxT1 ; . . . ; xTjIj Þ, uT ¼ ðuT1 ; . . . ; uTjIj Þ
454
J. Blomvall, P.O. Lindberg / European Journal of Operational Research 143 (2002) 452–461
and y T ¼ ðxT ; uT Þ. Furthermore, let Y P be the set of y that satisfy (1c)–(1e) and Y ¼ the set of y that satisfy constraint (1b). The T set of feasible solutions to (P) then is Y ¼ Y P Y ¼ . To allow for application of interior T point methods, we will assume that intðY P Þ Y ¼ is nonempty.
ðPl Þ min xi ;ui
def
fl ¼
X
½Ui ðxi ; ui Þ þ Ti ðl; xi ; ui Þ
ð2aÞ
i2I
s:t: xi ¼ Ai xi þ Bi ui þ bi ;
i 2 I n f1g;
ð2bÞ
where 2. The algorithm
Ti ðl; xi ; ui Þ ¼
2.1. Outline We will essentially apply a standard primal barrier method to problem (P). In the barrier subproblem for a given barrier penalty, Newton steps are determined by solving a second order expansion of the subproblem. In our case we solve this QP by viewing it as a LQCP over the scenario tree. The standard dynamic programming approach for discrete time LQCPs (e.g. Bertsekas, 1995, Chapter 4.1) carries over to the current situation. Thus we compute control laws that are linear in the states, by solving matrix Riccati equations backward in time. Applying these control laws, we compute new states, working forward in time. 2.2. Reformulation To initialize the interior point method we need a starting point that satisfies the inequality constraints strictly. In Appendix A we do this initialization. Some additional control variables are added and these give additional terms in the objective function, namely 1 2 ðIi ui gi Þ : 2l In the sequel we assume that the necessary modifications have been made to the original matrices and variables in (P). 2.3. Interior point formulation As is standard for primal barrier methods, the inequality constraints (1c)–(1e) are relaxed through a logarithmic barrier with parameter l. The barrier subproblem thus is
1 2 ðIi ui gi Þ 2l l1T lnðEi xi þ Fi ui ei Þ l1T lnðxi xli Þ l1T lnðxui xi Þ l1T lnðui uli Þ l1T lnðuui ui Þ:
Here we assume that ln operates componentwise on vectors and that 1 is the column unity vector of all ones. Further define def
fi ðxi ; ui Þ ¼ Ui ðxi ; ui Þ þ Ti ðl; xi ; ui Þ: 2.4. Second order approximation To solve (Pl ), we take Newton steps. To compute these, we make a second order approximation of (Pl ) at the current point xi ; ui . The gradients of the objective function, which are considered as row vectors, are def
T 1 Ei Eixi þ Fi ui ei !T !T 1 1 l þl u xi xli xi xi
qTi ¼ rxi fi ¼ rxi Ui l
and 1 def T riT ¼ rui fi ¼ rui Ui þ ðIi ui gi Þ Ii l T 1 l Fi Eixi þ Fi ui ei !T !T 1 1 l þl u : ui uli ui ui The vector divisions 1=ð Þ are made componentwise. The Hessians of the objective are
J. Blomvall, P.O. Lindberg / European Journal of Operational Research 143 (2002) 452–461 def
Qi ¼ r2x2 fi ¼ r2x2 Ui i
i
T þ lEi diag
2
1 Eixi þ Fi ui e i !!2 1 þ l diag xi xli !!2 1 þ l diag u ; xi xi
Ei
!
r2x2 Ui
r2xi ui Ui
r2ui xi Ui
r2u2 Ui þ l1 IiT Ii
i
455
;
i
is positive semidefinite. The second part, derived from the barrier term of (1c), is also positive semidefinite. The third, and final, part comes from the boundary constraints. This diagonal matrix has strictly positive diagonal elements and is hence positive definite. Thus Qi is positive definite. Corollary 2. Qi and Ri are positive definite.
def
Pi ¼ r2xi ui fi ¼ r2xi ui Ui þ lEiT diag
1 Eixi þ Fi ui e i
2
2.5. Newton step subproblem Fi ; When we have a solution xi ; ui (not necessarily feasible) to (Pl ) then xi ¼ Aixi þ Bi ui þ bi ð3Þ hold for appropriate bi . Using (3) we can now rewrite (2b)
T r2ui xi fi ¼ r2xi ui fi ¼ r2ui xi Ui þ lFiT diag
1 Eixi þ Fi ui e i
1 def Ri ¼ r2u2 fi ¼ r2u2 Ui þ IiT Ii i i l þ lFiT diag
2
1 Eixi þ Fi ui e i !!2 1 þ l diag ui uli !!2 1 þ l diag u : ui ui
Ei ;
Dxi ¼ Ai Dxi þ Bi Dui þ Dbi ; where Dbi ¼ bi bi . The second order expansion of (Pl ) is tree structured. It can be formulated as
2 Fi
The second order approximation of fi ðxi ; ui Þ at xi ; ui can be described as a function of the increment in the variables: def fi ðxi þ Dxi ; ui þ Dui Þ fi ðDxi ; Dui Þ ¼ ai þ qTi Dxi
þ 12DxTi Qi Dxi þ riT Dui þ 12DuTi Ri Dui þ DxTi Pi Dui ; where ai ¼ fi ðxi ; ui Þ. Proposition 1. Qi ¼
Qi PiT
Pi Ri
is positive definite.
ðPl Þ min
Dxi ;Dui
X
ai þ qTi Dxi þ 12DxTi Qi Dxi þ riT Dui
i2I
þ 12DuTi Ri Dui þ DxTi Pi Dui s:t: Dxi ¼ Ai Dxi þ Bi Dui þ Dbi ;
i 2 I n f1g:
ð4bÞ
This problem is tree structured in the sense that the controls at a given node i only influence states and constraints in succeeding nodes, j 2 Iþþ i . 2.5.1. Dynamic programming solution We solve (Pl ) by dynamic programming. To this end, we define value functions Vi ðDxi Þ ¼ optimal cost in ðPl Þ to terminal time from state Dxi at node i X aj þ qTj Dxj þ 12DxTj Qj Dxj ¼ min Duj ;j2J
j2J
þ
Proof. Qi can be divided into the sum of several parts. The first part,
ð4aÞ
X
rjT Duj þ 12DuTj Rj Duj þ DxTj Pj Duj
j2J
s:t: Dxj ¼ Aj Dxj þ Bj Duj þ Dbj ;
j 2 Iþþ i ;
456
J. Blomvall, P.O. Lindberg / European Journal of Operational Research 143 (2002) 452–461
S where J ¼ fig Iþþ i . As is standard in dynamic programming, the value functions satisfy the following recursion: Vi ðDxi Þ ¼ min fid ðDxi ; Dui ; fDxiþ giþ 2Iþ Þ Dui
¼ ai þ qTi Dxi þ 12DxTi Qi Dxi þ riT Dui
X
e i ¼ Qi þ Q Pei ¼ Pi þ
s:t: 8iþ 2 Iþ i :
ATiþ ðwiþ þ Wiþ Dbiþ Þ;
ð7dÞ
ð5bÞ
(Note that the dynamics constraint (5b) has been written in terms of Dxiþ instead of Dxi .) We make the following
X
ATiþ Wiþ Aiþ ;
ð7eÞ
iþ 2Iþ i
ð5aÞ
iþ 2Iþ i
Dxiþ ¼ Aiþ Dxi þ Biþ Dui þ Dbiþ
ð7cÞ
iþ 2Iþ i
def
þ 12DuTi Ri Dui þ DxTi Pi Dui X þ Viþ ðDxiþ Þ;
BTiþ Wiþ Biþ ;
iþ 2Iþ i
q~i ¼ qi þ
i
X
e i ¼ Ri þ R
X
ATiþ Wiþ Biþ :
ð7fÞ
iþ 2Iþ i
With these definitions the value function (5a) becomes Vi ðDxi Þ ¼ min fiu ðDxi ; Dui Þ Dui
def
e i Dxi ¼ a~i þ q~Ti Dxi þ 12DxTi Q T Te e i Dui : þ ð~ r þ Dx P i ÞDui þ 1DuT R
Inductive assumption. Vi is quadratic in Dxi , Vi ðDxi Þ ¼ ai þ wTi Dxi þ 12DxTi Wi Dxi ; where Wi is positive semidefinite.
Viþ ðDxiþ Þ ¼ aiþ þ
ðwTiþ
þ
þ
þ
1 DbTiþ Wiþ Dbiþ 2
DbTiþ Wiþ ÞAiþ Dxi
þ ðwTiþ þ DbTiþ Wiþ þ DxTi ATiþ Wiþ ÞBiþ Dui : To continue, the following notations are used: X a~i ¼ ai þ ðaiþ þ wTiþ Dbiþ þ 12DbTiþ Wiþ Dbiþ Þ;
Proof. By the inductive assumption Wiþ is positive semidefinite. Thus BTiþ Wiþ Biþ is positive semidefinite. Together with the fact that Ri is positive e i is positive definite. definite it follows that R By the lemma, the optimization problem in (8) is strictly convex, and the optimum can be found by differentiating fiu ðDxi ; Dui Þ with respect to Dui .
e i Dxi Vi ðDxi Þ ¼ a~i þ q~Ti Dxi þ 12DxTi Q e 1 ð~ri þ Pe T Dxi Þ ð~riT þ DxTi Pei Þ R i i T e 1 þ 12ð~ ri þ PeiT Dxi Þ R ri þ PeiT Dxi Þ i ð~
¼ ai þ wTi Dxi þ 12DxTi Wi Dxi ;
ð7aÞ
iþ 2Iþ i
ð9Þ
i
Inserting this optimal Dui in (8), Vi becomes
iþ 2Iþ i
BTiþ ðwiþ þ Wiþ Dbiþ Þ;
ð8Þ
i
Lemma 3. Under the inductive assumption, R~i is positive definite.
i
þ 12DuTi BTiþ Wiþ Biþ Dui
X
2
ei 0 :¼ rDui fiu ¼ ~riT þ DxTi Pei þ DuTi R e 1 ð~ri þ Pe T Dxi Þ: () D^ ui ¼ R
þ 12DxTi ATiþ Wiþ Aiþ Dxi
~ ri ¼ ri þ
i
2.5.2. Optimal control laws
The inductive assumption holds true for all successors of i 2 T, since in these nodes ai ¼ 0, wi ¼ 0 and Wi ¼ 0. Hence we assume it true for all immediate successors to an arbitrary node i, and shall prove it true for i. Thereby proving it for all nodes. In the process we will prove that the controls Dui are affine in the state Dxi and derive the recursions that give these relations. Using the expression (5b) for Dxiþ we can rewrite the value function (6) for a successor node. wTiþ Dbiþ
i
ð6Þ
where ð7bÞ
e 1 r~i ; ai ¼ a~i 12~riT R i
ð10aÞ
J. Blomvall, P.O. Lindberg / European Journal of Operational Research 143 (2002) 452–461
e 1 ~ wi ¼ q~i Pei R i ri ;
ð10bÞ
e i Pei R e 1 Pe T : Wi ¼ Q i i
ð10cÞ
457
cessors of i 2 T. Thus it holds for all i 2 I by induction. Note that we now also have the equations to determine the matrices in (6).
The recursion (10c) (together with Eqs. (7a)– (7f)) is known as the (discrete time) Riccati equation in control theory.
2.5.3. Summing up We now have the complete algorithm to solve (Pl ).
Lemma 4. If f ðy; zÞ is a strictly convex quadratic function in y, convex quadratic in z and separable in y, z, then f ðy; Ay þ bÞ is strictly convex quadratic.
Theorem 7. The Newton step determination problem ðPl Þ for the barrier subproblem ðPl Þ is solved by (9), where the matrices, vectors and constants are defined recursively by (7a)–(7f) and (10a)–(10c).
Proof. The sum of one strictly convex function and one convex function is a strictly convex function. That f ðy; Ay þ bÞ is quadratic is immediate.
Proof. This follows from the previous analysis.
Proof. Strict convexity follows since f ðy; Ay þ bÞ is a restrification of f ðy; zÞ to a subspace. That f ðy; Ay þ bÞ is quadratic is immediate.
The solution procedure of (Pl ) now works as follows. Determine control laws that are linear in the state by calculating (7a)–(7f) and (10a)–(10c) recursively starting at the terminal nodes. When the root node is reached, we can then using the initial state x1 , the control laws (9) and the dynamics constraint (4b) determine the optimal solution, by making a sweep from the root node to the terminal nodes.
Theorem 6. Wi is a positive definite matrix if Wiþ is positive semidefinite for iþ 2 Iþ i .
2.6. Solution procedure
Lemma 5. If f ðy; zÞ is a strictly convex quadratic function, then f ðy; Ay þ bÞ is strictly convex quadratic.
Proof. The function fid ðDxi ; Dui ; fDxiþ giþ 2Iþ Þ dei fined in (5a) is convex quadratic being the sum of convex quadratic functions of the variables in each node in question. From the definition of fiu ðDxi ; Dui Þ (8) we see that fiu ðDxi ; Dui Þ ¼ fid ðDxi ; Dui ; fAiþ Dxi þ Biþ Dui þ Dbiþ giþ 2Iþ Þ: i
Applying Lemma 4 we conclude that fiu ðDxi ; Dui Þ is strictly convex quadratic. Finally we have that uiÞ Vi ðDxi Þ ¼ fiu ðDxi ; Db e 1 ½~ eT ¼ fiu ðDxi ; R i ri þ P i Dxi Þ: Using Lemma 5 it follows that Vi ðDxi Þ is strictly convex quadratic, thus Wi is positive definite. By Theorem 6, if (6) holds for all iþ 2 Iþ i then it also holds for i. As noted (6) holds for all suc-
In iteration j a second order approximation is made at the point yj 2 intðY P Þ, for a parameter lj . Theorem 7 implies that one can determine the optimal solution yj þ Dyj to (Pl ), by first determining optimal control laws and then, since we know x1 , the optimal decisions and states. We know that yj þ Dyj 2 Y ¼ . We might however violate the inequality constraints. If aj is the largest a such that yj þ aDyj 2 Y P , then we select the step length as aj ¼ minðnaj ; 1Þ where n 2 ð0; 1Þ. This gives the modified Newton step. With this step length the new solution is yjþ1 ¼ yj þ aj Dyj 2 Y P . l is decreased each iteration, until we reach a target level l. Then l is kept fixed and we iterate until there is no significant improvement in the barrier subproblem objective (after having had a quadratic convergence rate). Then, we assume that we are sufficiently close to the central path for a sufficiently small l. The solution procedure is illustrated in Fig. 2.
458
J. Blomvall, P.O. Lindberg / European Journal of Operational Research 143 (2002) 452–461
amount xai , owned in the previous node, i . Furthermore we know the amounts bought, uab i , and sold, uas i . The total amount in the current node thus is as xai ¼ xai þ uab i ui
a 2 Ai ;
where as xai ; uab i ; ui P 0:
3.1.2. Cash Similarly to when the asset constraints were defined, we know that in the previous node each asset was bought at the price piab and sold at pias . If the cash held in the previous node was xli , then the following holds for the cash in the current node: ! X l l as as ab ab xi ¼ xi þ ðpi ui pi ui Þ erl ; a2Ai
xli
Fig. 2. Solution method.
3. Computational results We apply the method above to a financial application. It is a multistage stochastic program where we maximize the expected utility of a call option portfolio. We use one underlying asset, the OMX-index, a Swedish composite index similar to the S&P 500 index. The other assets are call options on the OMX-index and one risk-free asset. This is an extension of the two stage model used in Blomvall and Lindberg (in press). 3.1. Model formulation By making a mathematical formulation of the investment problem, assuming discrete prices we get a formulation that fits (P). 3.1.1. Assets In node i the available assets are represented by the set Ai . Available assets can differ between the nodes (since options have a limited lifespan). For each node i and each asset a 2 Ai we know the
where P 0 is required and rl is the one period lending rate (we assume a fixed length of the periods). 3.1.3. Objective function Utility theory has been treated extensively in economic theory. It is possible to show, under some few reasonable assumptions on the preferences of a rational investor, that he/she will act so as to maximize some expected utility (e.g. von Neumann and Morgenstern, 1953). The objective function, X max pi U ðcTi xi Þ;
ui ;i2I
i2T
corresponds to utility maximization if cTi xi denotes the total value of the portfolio in node i and pi denotes the probability of ending up in that node. According to Kelly (1956) there exists one special utility function (the logarithmic utility) which outperforms any other utility function in the long run in terms of growth of wealth (under some natural assumptions). If an investor invests using the logarithmic utility function then he/she will with probability 1 have more money in the long run than if any other utility function had been used. Therefore we will use the logarithmic utility
J. Blomvall, P.O. Lindberg / European Journal of Operational Research 143 (2002) 452–461
459
function in our application. This approach has given very good result in our previous application (Blomvall and Lindberg, in press). From the above, it is clear that a multistage asset investment problem with discrete outcomes of prices can be modeled by (P). 3.2. Settings The same method as in (Blomvall and Lindberg, in press) is used to estimate the probabilites of the scenarios. The option prices in the nodes are estimated using so called implied volatilities. In the root node we use market prices. We initiate the model on February 23, 1990. On this day there were 10 call options that matured during the next month. The transaction costs for trading in the OMX-index and call options are 0.1% and 2.0%, respectively. We use a seven stage model, with 10 outcomes in each node. This yields 1,000,000 scenarios and 1,111,111 nodes in the tree. The resulting optimization problem has 15,777,740 variables and 16,888,850 constraints. In the solver we used the parameter values n ¼ 0:95 and l ¼ 1016 . To solve the problem we use a 200 MHz Sun Sparc station with 2 GB memory. It takes approximately 12 minutes to solve one subproblem (Pl ). The overall problem is solved in approximately 100 iterations (20 hours), depending on the termination criterion. As a comparison Gondzio and Kouwenberg (2001) solve a linear optimization problem with 4,826,809 scenarios, 24,938,502 variables and 12,469,250 constraints in 5 hours, but on a Parsitech CC16 with 13 processors. As can be seen in Fig. 3, the objective function value of (P) decreases (max problem) in the first 20 iterations. Since we only take one Newton step, before changing l, it takes 20 iterations to reach the close proximity of the central path. When the central path has been reached the objective function value is steadily increasing. In Fig. 4 we are instead looking at the rate of convergence towards the ‘‘optimal’’ objective function value. We use the objective function value of (P) in iteration 192 as an approximation of the optimal function value. One interesting observation is that the deviation
Fig. 3. Convergence of the objective function (maximization problem).
from the optimal value decreases with a constant factor in each iteration.
4. Summary and outlook The dual block angular structure in stochastic programming problems has for a long time been used to decompose the optimization problem. Interior point methods has been used in recent years as a means to solve the problem, either directly (Birge and Holmes, 1992) or as a solver in the subproblems (Gondzio and Kouwenberg, 2001). Our approach is a standard primal barrier approach. Since the stochastic programming problem was stated in control form with bounds, we get a very convenient direction finding subproblem: An optimization problem with strictly convex quadratic objective function and linear equality constraints. This optimization problem is solved by an equation system; we only need an efficient way to do this. As is standard practice in control theory (Bertsekas, 1995), the subproblem is solved by dynamic programming. The control theory approach, based on Riccati equations, carries over to the present tree structured situation (see also Salinger and Rockafellar, 1999). For the numerical results, we used an extended version of the financial model in (Blomvall and Lindberg, in press). With the solver we managed to solve a problem with
460
J. Blomvall, P.O. Lindberg / European Journal of Operational Research 143 (2002) 452–461
Fig. 4. Convergence of the objective function towards the optimal objective funtion value.
nonlinear separable objective function, 15,777,740 variables and 16,888,850 constraints in 20 hours. The number of scenarios was 1,000,000. Thus we have developed a solver that can solve very large scale multistage stochastic programming problems that have linear or nonlinear separable objective function. The approach presented is quite versatile. In a forthcoming paper, we will give an extension to more general stochastic programming problems. Further, the approach obviously is easily parallelizable. We are in the process to develop a parallel version for a PC-cluster.
We now have to modify the matrices that premultiply the control variables. For the constraint (1b) we make the following modifications:
Appendix A. Initialization
where M is a large positive value. In order to find a starting solution we first initialize ui in a way to make as many as possible of the constraints (1c) feasible. It is only necessary to introduce new control variables, uiP , for those constraints that are not satisfied. This reduces the number of new variables. We want the recently introduced variables, uiP , to become negative, so that they can be eliminated from the optimization problem. To achieve negative uiP the following term is added to the objective function:
We have to initialize the algorithm if we do not have a feasible solution to start from. We want to find a y that satisfies (1c)–(1e) in order to initialize the interior point method. To find a solution that satisfies (1c) we introduce additional control variables, uiP . To simplify the notation uiP has the same dimension as ei . This gives us a new vector of control variables ðu0i ÞT ¼ uTi ; ðuiP ÞT i 2 I:
B0i ¼ ðBi ; 0Þ i 2 I n f1g: In the inequality constraint (1c), we use the variable uiP to find a feasible solution. This leads to the following modifications: Fi0 ¼ ðFi ; I Þ
i 2 I;
where I denotes the identity matrix. Finally (1e) is modified as
0 uli ¼ uli ; 1M i 2 I;
0 uui ¼ uui ; 1M i 2 I;
J. Blomvall, P.O. Lindberg / European Journal of Operational Research 143 (2002) 452–461 1 ðI u0 2l i i
2
gi Þ ;
where 0 0 Ii ¼ 0 I
i2I
gi is the target values for uiP and l is a parameter. References Bertsekas, D.P., 1995. Dynamic Programming and Optimal Control. Athena Scientific, Belmont, MA. Birge, J.R., 1997. Stochastic programming computation and applications. INFORMS Journal on Computing 9 (2), 111– 133. Birge, J., Holmes, D., 1992. Efficient solution of two-stage stochastic linear programs using interior point methods. Computational Optimization and Applications 1 (3), 245–276. Blomvall, J., Lindberg, P.O. Back-testing the performance of an actively managed option portfolio at the Swedish stock
461
market, 1990-1999, Journal of Economic Dynamics and Control, in press. Gondzio, J., Kouwenberg, R., 2001. High performance computing for asset liability management. Operations Research 49, 879–891. Kelly, J.L., 1956. A new interpretation of information rate. Bell System Technical Journal 35, 917–926. Salinger, D.H., Rockafellar, R.T., 1999. Dynamic splitting: An algorithm for deterministic and stochastic multiperiod optimization, working paper, Department of Mathematics, University of Washington, Seattle. Steinbach, M.C., 2000. Hierarchical sparsity in multistage convex stochastic programs, Report 00-15, ZIB, Berlin. Van Slyke, R., Wets, R., 1969. L-shaped linear programs with applications to control and stochastic programming. SIAM Journal on Applied Mathematics 17, 638– 663. Vladimirou, H., Zenios, S.A., 1999. Scalable parallel computations for large-scale stochastic programming. Annals of Operations Research 90, 87–129. von Neumann, J., Morgenstern, O., 1953. Theory of Games and Economic Behavior. Princeton University Press, Princeton.
A Riccati-based primal interior point solver for multistage stochastic programming €rgen Blomvall *, Per Olov Lindberg Jo Department of Mathematics, Link€opings Universitet, SE-58183 Link€oping, Sweden Received 11 January 2000; accepted 12 July 2001
Abstract We propose a new method for certain multistage stochastic programs with linear or nonlinear objective function, combining a primal interior point approach with a linear-quadratic control problem over the scenario tree. The latter problem, which is the direction finding problem for the barrier subproblem is solved through dynamic programming using Riccati equations. In this way we combine the low iteration count of interior point methods with an efficient solver for the subproblems. The computational results are promising. We have solved a financial problem with 1,000,000 scenarios, 15,777,740 variables and 16,888,850 constraints in 20 hours on a moderate computer. Ó 2002 Elsevier Science B.V. All rights reserved. Keywords: Dynamic programming; Finance; Interior point methods; Stochastic programming
1. Introduction Stochastic programming has many important applications (see e.g. Birge, 1997). It further easily leads to huge mathematical programs. It would thus be a natural candidate for the application of interior point methods. However, a direct application typically leads to dense systems of normal equations, essentially preventing this approach. Some authors have worked on circumventing these difficulties (Birge and Holmes, 1992) including work aiming at parallel computations (Vladimirou and Zenios, 1999). See also Birge (1997) and the references therein. *
Corresponding author. Tel.: +46-13-281406; fax: +46-13285770. E-mail address: [email protected] (J. Blomvall).
Others (e.g. Gondzio and Kouwenberg, 2001) have used interior point methods to solve the problems that appear in decomposition schemes for stochastic programming, such as L-shaped decompostion (Van Slyke and Wets, 1969). In the current paper we apply a primal interior point method to certain stochastic programs. Our main source of motivation comes from the financial area. Our work is related to, but independent of that of Steinbach (2000), whose approach is based on primal–dual interior point methods. Our approach is based on viewing the problem of computing a Newton step in the barrier subproblem as a linear-quadratic control problem (LQCP) over the scenario tree. For standard LQCPs there are standard techniques based on Riccati equations to compute closed form solutions (see e.g. Bertsekas, 1995). These techniques carry over
0377-2217/02/$ - see front matter Ó 2002 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 7 - 2 2 1 7 ( 0 2 ) 0 0 3 0 1 - 6
J. Blomvall, P.O. Lindberg / European Journal of Operational Research 143 (2002) 452–461
to our situation as we will show. Salinger and Rockafellar (1999) in applying operator splitting to stochastic programs also arrive at LQCPs over trees and observe that they can be treated by Riccati equations. In the rest of Section 1, we will give our model framework and state our model. Section 2 is devoted to deriving the algorithm. Our test problem, which motivated this study, is given in Section 3 together with some computational results. We are able to solve quite large problems (1,000,000 scenarios giving 15,777,740 nontrivial variables and 16,888,850 constraints) in reasonable time (20 hours) with rather limited computer resources. 1.1. Model framework Stochastic programs with discrete probability distribution can equivalently be stated as mathematical programs. These can be represented by time-staged trees, scenario trees (see Fig. 1), where each node represents a certain outcome of the random variables up to that stage. The stages (the depth from the root node) represent the evolution of time. Each node i 2 I has one predecessor
node, i , and a set, Iþ i , of immediate successor nodes in the next stage. (This is the same notation as in Salinger and Rockafellar (1999).) The leftmost node, f1g, in the tree is the root node. In the first time stage there is no uncertainty. The terminal nodes, T, represent the last time stage. The set Iþþ contains all the successors of node i, i.e. all i successors, up to and including the leaf nodes. 1.2. Model The notation above allows for convenient formulation of many stochastic programming problems, for example in finance which has been our main source of inspiration. Assume that we want to solve the following stochastic programming problem (P), partly in ‘‘control form’’: ðP Þ min xi ;ui
X
Ui ðxi ; ui Þ
Ei x i þ F i u i P e i ; ui 2
ð1aÞ
i2I
s:t: xi ¼ Ai xi þ Bi ui þ bi ; xi 2
Fig. 1. Notation for a multistage model.
453
½xli ; xui ; ½uli ; uui ;
i 2 I;
i 2 I n f1g;
ð1bÞ ð1cÞ
i 2 I;
ð1dÞ
i 2 I;
ð1eÞ
where x1 is a constant. We say that the problem is in control form, since the decision variables, are partitioned into state variables xi and control variables ui . For any values of the controls fui gi2I the states fxi gi2I are uniquely determined through state transitions equations (1b) (note that (1b) is of the same format as in a discrete time control problem). On top of (1b), we have local (node) constraints (1c) and simple bounds (1d)–(1e). The objective is a sum of convex functions of the node variables. The probability of each node is assumed to be included in Ui ð Þ. Many stochastic programs can be formulated as (P). In our main financial application, the ui are buying and selling decisions, whereas the xi are ‘‘inventories’’ of stocks, options, cash etc. Letting j j denote the number of elements of a set then we can define the following composite variables: xT ¼ ðxT1 ; . . . ; xTjIj Þ, uT ¼ ðuT1 ; . . . ; uTjIj Þ
454
J. Blomvall, P.O. Lindberg / European Journal of Operational Research 143 (2002) 452–461
and y T ¼ ðxT ; uT Þ. Furthermore, let Y P be the set of y that satisfy (1c)–(1e) and Y ¼ the set of y that satisfy constraint (1b). The T set of feasible solutions to (P) then is Y ¼ Y P Y ¼ . To allow for application of interior T point methods, we will assume that intðY P Þ Y ¼ is nonempty.
ðPl Þ min xi ;ui
def
fl ¼
X
½Ui ðxi ; ui Þ þ Ti ðl; xi ; ui Þ
ð2aÞ
i2I
s:t: xi ¼ Ai xi þ Bi ui þ bi ;
i 2 I n f1g;
ð2bÞ
where 2. The algorithm
Ti ðl; xi ; ui Þ ¼
2.1. Outline We will essentially apply a standard primal barrier method to problem (P). In the barrier subproblem for a given barrier penalty, Newton steps are determined by solving a second order expansion of the subproblem. In our case we solve this QP by viewing it as a LQCP over the scenario tree. The standard dynamic programming approach for discrete time LQCPs (e.g. Bertsekas, 1995, Chapter 4.1) carries over to the current situation. Thus we compute control laws that are linear in the states, by solving matrix Riccati equations backward in time. Applying these control laws, we compute new states, working forward in time. 2.2. Reformulation To initialize the interior point method we need a starting point that satisfies the inequality constraints strictly. In Appendix A we do this initialization. Some additional control variables are added and these give additional terms in the objective function, namely 1 2 ðIi ui gi Þ : 2l In the sequel we assume that the necessary modifications have been made to the original matrices and variables in (P). 2.3. Interior point formulation As is standard for primal barrier methods, the inequality constraints (1c)–(1e) are relaxed through a logarithmic barrier with parameter l. The barrier subproblem thus is
1 2 ðIi ui gi Þ 2l l1T lnðEi xi þ Fi ui ei Þ l1T lnðxi xli Þ l1T lnðxui xi Þ l1T lnðui uli Þ l1T lnðuui ui Þ:
Here we assume that ln operates componentwise on vectors and that 1 is the column unity vector of all ones. Further define def
fi ðxi ; ui Þ ¼ Ui ðxi ; ui Þ þ Ti ðl; xi ; ui Þ: 2.4. Second order approximation To solve (Pl ), we take Newton steps. To compute these, we make a second order approximation of (Pl ) at the current point xi ; ui . The gradients of the objective function, which are considered as row vectors, are def
T 1 Ei Eixi þ Fi ui ei !T !T 1 1 l þl u xi xli xi xi
qTi ¼ rxi fi ¼ rxi Ui l
and 1 def T riT ¼ rui fi ¼ rui Ui þ ðIi ui gi Þ Ii l T 1 l Fi Eixi þ Fi ui ei !T !T 1 1 l þl u : ui uli ui ui The vector divisions 1=ð Þ are made componentwise. The Hessians of the objective are
J. Blomvall, P.O. Lindberg / European Journal of Operational Research 143 (2002) 452–461 def
Qi ¼ r2x2 fi ¼ r2x2 Ui i
i
T þ lEi diag
2
1 Eixi þ Fi ui e i !!2 1 þ l diag xi xli !!2 1 þ l diag u ; xi xi
Ei
!
r2x2 Ui
r2xi ui Ui
r2ui xi Ui
r2u2 Ui þ l1 IiT Ii
i
455
;
i
is positive semidefinite. The second part, derived from the barrier term of (1c), is also positive semidefinite. The third, and final, part comes from the boundary constraints. This diagonal matrix has strictly positive diagonal elements and is hence positive definite. Thus Qi is positive definite. Corollary 2. Qi and Ri are positive definite.
def
Pi ¼ r2xi ui fi ¼ r2xi ui Ui þ lEiT diag
1 Eixi þ Fi ui e i
2
2.5. Newton step subproblem Fi ; When we have a solution xi ; ui (not necessarily feasible) to (Pl ) then xi ¼ Aixi þ Bi ui þ bi ð3Þ hold for appropriate bi . Using (3) we can now rewrite (2b)
T r2ui xi fi ¼ r2xi ui fi ¼ r2ui xi Ui þ lFiT diag
1 Eixi þ Fi ui e i
1 def Ri ¼ r2u2 fi ¼ r2u2 Ui þ IiT Ii i i l þ lFiT diag
2
1 Eixi þ Fi ui e i !!2 1 þ l diag ui uli !!2 1 þ l diag u : ui ui
Ei ;
Dxi ¼ Ai Dxi þ Bi Dui þ Dbi ; where Dbi ¼ bi bi . The second order expansion of (Pl ) is tree structured. It can be formulated as
2 Fi
The second order approximation of fi ðxi ; ui Þ at xi ; ui can be described as a function of the increment in the variables: def fi ðxi þ Dxi ; ui þ Dui Þ fi ðDxi ; Dui Þ ¼ ai þ qTi Dxi
þ 12DxTi Qi Dxi þ riT Dui þ 12DuTi Ri Dui þ DxTi Pi Dui ; where ai ¼ fi ðxi ; ui Þ. Proposition 1. Qi ¼
Qi PiT
Pi Ri
is positive definite.
ðPl Þ min
Dxi ;Dui
X
ai þ qTi Dxi þ 12DxTi Qi Dxi þ riT Dui
i2I
þ 12DuTi Ri Dui þ DxTi Pi Dui s:t: Dxi ¼ Ai Dxi þ Bi Dui þ Dbi ;
i 2 I n f1g:
ð4bÞ
This problem is tree structured in the sense that the controls at a given node i only influence states and constraints in succeeding nodes, j 2 Iþþ i . 2.5.1. Dynamic programming solution We solve (Pl ) by dynamic programming. To this end, we define value functions Vi ðDxi Þ ¼ optimal cost in ðPl Þ to terminal time from state Dxi at node i X aj þ qTj Dxj þ 12DxTj Qj Dxj ¼ min Duj ;j2J
j2J
þ
Proof. Qi can be divided into the sum of several parts. The first part,
ð4aÞ
X
rjT Duj þ 12DuTj Rj Duj þ DxTj Pj Duj
j2J
s:t: Dxj ¼ Aj Dxj þ Bj Duj þ Dbj ;
j 2 Iþþ i ;
456
J. Blomvall, P.O. Lindberg / European Journal of Operational Research 143 (2002) 452–461
S where J ¼ fig Iþþ i . As is standard in dynamic programming, the value functions satisfy the following recursion: Vi ðDxi Þ ¼ min fid ðDxi ; Dui ; fDxiþ giþ 2Iþ Þ Dui
¼ ai þ qTi Dxi þ 12DxTi Qi Dxi þ riT Dui
X
e i ¼ Qi þ Q Pei ¼ Pi þ
s:t: 8iþ 2 Iþ i :
ATiþ ðwiþ þ Wiþ Dbiþ Þ;
ð7dÞ
ð5bÞ
(Note that the dynamics constraint (5b) has been written in terms of Dxiþ instead of Dxi .) We make the following
X
ATiþ Wiþ Aiþ ;
ð7eÞ
iþ 2Iþ i
ð5aÞ
iþ 2Iþ i
Dxiþ ¼ Aiþ Dxi þ Biþ Dui þ Dbiþ
ð7cÞ
iþ 2Iþ i
def
þ 12DuTi Ri Dui þ DxTi Pi Dui X þ Viþ ðDxiþ Þ;
BTiþ Wiþ Biþ ;
iþ 2Iþ i
q~i ¼ qi þ
i
X
e i ¼ Ri þ R
X
ATiþ Wiþ Biþ :
ð7fÞ
iþ 2Iþ i
With these definitions the value function (5a) becomes Vi ðDxi Þ ¼ min fiu ðDxi ; Dui Þ Dui
def
e i Dxi ¼ a~i þ q~Ti Dxi þ 12DxTi Q T Te e i Dui : þ ð~ r þ Dx P i ÞDui þ 1DuT R
Inductive assumption. Vi is quadratic in Dxi , Vi ðDxi Þ ¼ ai þ wTi Dxi þ 12DxTi Wi Dxi ; where Wi is positive semidefinite.
Viþ ðDxiþ Þ ¼ aiþ þ
ðwTiþ
þ
þ
þ
1 DbTiþ Wiþ Dbiþ 2
DbTiþ Wiþ ÞAiþ Dxi
þ ðwTiþ þ DbTiþ Wiþ þ DxTi ATiþ Wiþ ÞBiþ Dui : To continue, the following notations are used: X a~i ¼ ai þ ðaiþ þ wTiþ Dbiþ þ 12DbTiþ Wiþ Dbiþ Þ;
Proof. By the inductive assumption Wiþ is positive semidefinite. Thus BTiþ Wiþ Biþ is positive semidefinite. Together with the fact that Ri is positive e i is positive definite. definite it follows that R By the lemma, the optimization problem in (8) is strictly convex, and the optimum can be found by differentiating fiu ðDxi ; Dui Þ with respect to Dui .
e i Dxi Vi ðDxi Þ ¼ a~i þ q~Ti Dxi þ 12DxTi Q e 1 ð~ri þ Pe T Dxi Þ ð~riT þ DxTi Pei Þ R i i T e 1 þ 12ð~ ri þ PeiT Dxi Þ R ri þ PeiT Dxi Þ i ð~
¼ ai þ wTi Dxi þ 12DxTi Wi Dxi ;
ð7aÞ
iþ 2Iþ i
ð9Þ
i
Inserting this optimal Dui in (8), Vi becomes
iþ 2Iþ i
BTiþ ðwiþ þ Wiþ Dbiþ Þ;
ð8Þ
i
Lemma 3. Under the inductive assumption, R~i is positive definite.
i
þ 12DuTi BTiþ Wiþ Biþ Dui
X
2
ei 0 :¼ rDui fiu ¼ ~riT þ DxTi Pei þ DuTi R e 1 ð~ri þ Pe T Dxi Þ: () D^ ui ¼ R
þ 12DxTi ATiþ Wiþ Aiþ Dxi
~ ri ¼ ri þ
i
2.5.2. Optimal control laws
The inductive assumption holds true for all successors of i 2 T, since in these nodes ai ¼ 0, wi ¼ 0 and Wi ¼ 0. Hence we assume it true for all immediate successors to an arbitrary node i, and shall prove it true for i. Thereby proving it for all nodes. In the process we will prove that the controls Dui are affine in the state Dxi and derive the recursions that give these relations. Using the expression (5b) for Dxiþ we can rewrite the value function (6) for a successor node. wTiþ Dbiþ
i
ð6Þ
where ð7bÞ
e 1 r~i ; ai ¼ a~i 12~riT R i
ð10aÞ
J. Blomvall, P.O. Lindberg / European Journal of Operational Research 143 (2002) 452–461
e 1 ~ wi ¼ q~i Pei R i ri ;
ð10bÞ
e i Pei R e 1 Pe T : Wi ¼ Q i i
ð10cÞ
457
cessors of i 2 T. Thus it holds for all i 2 I by induction. Note that we now also have the equations to determine the matrices in (6).
The recursion (10c) (together with Eqs. (7a)– (7f)) is known as the (discrete time) Riccati equation in control theory.
2.5.3. Summing up We now have the complete algorithm to solve (Pl ).
Lemma 4. If f ðy; zÞ is a strictly convex quadratic function in y, convex quadratic in z and separable in y, z, then f ðy; Ay þ bÞ is strictly convex quadratic.
Theorem 7. The Newton step determination problem ðPl Þ for the barrier subproblem ðPl Þ is solved by (9), where the matrices, vectors and constants are defined recursively by (7a)–(7f) and (10a)–(10c).
Proof. The sum of one strictly convex function and one convex function is a strictly convex function. That f ðy; Ay þ bÞ is quadratic is immediate.
Proof. This follows from the previous analysis.
Proof. Strict convexity follows since f ðy; Ay þ bÞ is a restrification of f ðy; zÞ to a subspace. That f ðy; Ay þ bÞ is quadratic is immediate.
The solution procedure of (Pl ) now works as follows. Determine control laws that are linear in the state by calculating (7a)–(7f) and (10a)–(10c) recursively starting at the terminal nodes. When the root node is reached, we can then using the initial state x1 , the control laws (9) and the dynamics constraint (4b) determine the optimal solution, by making a sweep from the root node to the terminal nodes.
Theorem 6. Wi is a positive definite matrix if Wiþ is positive semidefinite for iþ 2 Iþ i .
2.6. Solution procedure
Lemma 5. If f ðy; zÞ is a strictly convex quadratic function, then f ðy; Ay þ bÞ is strictly convex quadratic.
Proof. The function fid ðDxi ; Dui ; fDxiþ giþ 2Iþ Þ dei fined in (5a) is convex quadratic being the sum of convex quadratic functions of the variables in each node in question. From the definition of fiu ðDxi ; Dui Þ (8) we see that fiu ðDxi ; Dui Þ ¼ fid ðDxi ; Dui ; fAiþ Dxi þ Biþ Dui þ Dbiþ giþ 2Iþ Þ: i
Applying Lemma 4 we conclude that fiu ðDxi ; Dui Þ is strictly convex quadratic. Finally we have that uiÞ Vi ðDxi Þ ¼ fiu ðDxi ; Db e 1 ½~ eT ¼ fiu ðDxi ; R i ri þ P i Dxi Þ: Using Lemma 5 it follows that Vi ðDxi Þ is strictly convex quadratic, thus Wi is positive definite. By Theorem 6, if (6) holds for all iþ 2 Iþ i then it also holds for i. As noted (6) holds for all suc-
In iteration j a second order approximation is made at the point yj 2 intðY P Þ, for a parameter lj . Theorem 7 implies that one can determine the optimal solution yj þ Dyj to (Pl ), by first determining optimal control laws and then, since we know x1 , the optimal decisions and states. We know that yj þ Dyj 2 Y ¼ . We might however violate the inequality constraints. If aj is the largest a such that yj þ aDyj 2 Y P , then we select the step length as aj ¼ minðnaj ; 1Þ where n 2 ð0; 1Þ. This gives the modified Newton step. With this step length the new solution is yjþ1 ¼ yj þ aj Dyj 2 Y P . l is decreased each iteration, until we reach a target level l. Then l is kept fixed and we iterate until there is no significant improvement in the barrier subproblem objective (after having had a quadratic convergence rate). Then, we assume that we are sufficiently close to the central path for a sufficiently small l. The solution procedure is illustrated in Fig. 2.
458
J. Blomvall, P.O. Lindberg / European Journal of Operational Research 143 (2002) 452–461
amount xai , owned in the previous node, i . Furthermore we know the amounts bought, uab i , and sold, uas i . The total amount in the current node thus is as xai ¼ xai þ uab i ui
a 2 Ai ;
where as xai ; uab i ; ui P 0:
3.1.2. Cash Similarly to when the asset constraints were defined, we know that in the previous node each asset was bought at the price piab and sold at pias . If the cash held in the previous node was xli , then the following holds for the cash in the current node: ! X l l as as ab ab xi ¼ xi þ ðpi ui pi ui Þ erl ; a2Ai
xli
Fig. 2. Solution method.
3. Computational results We apply the method above to a financial application. It is a multistage stochastic program where we maximize the expected utility of a call option portfolio. We use one underlying asset, the OMX-index, a Swedish composite index similar to the S&P 500 index. The other assets are call options on the OMX-index and one risk-free asset. This is an extension of the two stage model used in Blomvall and Lindberg (in press). 3.1. Model formulation By making a mathematical formulation of the investment problem, assuming discrete prices we get a formulation that fits (P). 3.1.1. Assets In node i the available assets are represented by the set Ai . Available assets can differ between the nodes (since options have a limited lifespan). For each node i and each asset a 2 Ai we know the
where P 0 is required and rl is the one period lending rate (we assume a fixed length of the periods). 3.1.3. Objective function Utility theory has been treated extensively in economic theory. It is possible to show, under some few reasonable assumptions on the preferences of a rational investor, that he/she will act so as to maximize some expected utility (e.g. von Neumann and Morgenstern, 1953). The objective function, X max pi U ðcTi xi Þ;
ui ;i2I
i2T
corresponds to utility maximization if cTi xi denotes the total value of the portfolio in node i and pi denotes the probability of ending up in that node. According to Kelly (1956) there exists one special utility function (the logarithmic utility) which outperforms any other utility function in the long run in terms of growth of wealth (under some natural assumptions). If an investor invests using the logarithmic utility function then he/she will with probability 1 have more money in the long run than if any other utility function had been used. Therefore we will use the logarithmic utility
J. Blomvall, P.O. Lindberg / European Journal of Operational Research 143 (2002) 452–461
459
function in our application. This approach has given very good result in our previous application (Blomvall and Lindberg, in press). From the above, it is clear that a multistage asset investment problem with discrete outcomes of prices can be modeled by (P). 3.2. Settings The same method as in (Blomvall and Lindberg, in press) is used to estimate the probabilites of the scenarios. The option prices in the nodes are estimated using so called implied volatilities. In the root node we use market prices. We initiate the model on February 23, 1990. On this day there were 10 call options that matured during the next month. The transaction costs for trading in the OMX-index and call options are 0.1% and 2.0%, respectively. We use a seven stage model, with 10 outcomes in each node. This yields 1,000,000 scenarios and 1,111,111 nodes in the tree. The resulting optimization problem has 15,777,740 variables and 16,888,850 constraints. In the solver we used the parameter values n ¼ 0:95 and l ¼ 1016 . To solve the problem we use a 200 MHz Sun Sparc station with 2 GB memory. It takes approximately 12 minutes to solve one subproblem (Pl ). The overall problem is solved in approximately 100 iterations (20 hours), depending on the termination criterion. As a comparison Gondzio and Kouwenberg (2001) solve a linear optimization problem with 4,826,809 scenarios, 24,938,502 variables and 12,469,250 constraints in 5 hours, but on a Parsitech CC16 with 13 processors. As can be seen in Fig. 3, the objective function value of (P) decreases (max problem) in the first 20 iterations. Since we only take one Newton step, before changing l, it takes 20 iterations to reach the close proximity of the central path. When the central path has been reached the objective function value is steadily increasing. In Fig. 4 we are instead looking at the rate of convergence towards the ‘‘optimal’’ objective function value. We use the objective function value of (P) in iteration 192 as an approximation of the optimal function value. One interesting observation is that the deviation
Fig. 3. Convergence of the objective function (maximization problem).
from the optimal value decreases with a constant factor in each iteration.
4. Summary and outlook The dual block angular structure in stochastic programming problems has for a long time been used to decompose the optimization problem. Interior point methods has been used in recent years as a means to solve the problem, either directly (Birge and Holmes, 1992) or as a solver in the subproblems (Gondzio and Kouwenberg, 2001). Our approach is a standard primal barrier approach. Since the stochastic programming problem was stated in control form with bounds, we get a very convenient direction finding subproblem: An optimization problem with strictly convex quadratic objective function and linear equality constraints. This optimization problem is solved by an equation system; we only need an efficient way to do this. As is standard practice in control theory (Bertsekas, 1995), the subproblem is solved by dynamic programming. The control theory approach, based on Riccati equations, carries over to the present tree structured situation (see also Salinger and Rockafellar, 1999). For the numerical results, we used an extended version of the financial model in (Blomvall and Lindberg, in press). With the solver we managed to solve a problem with
460
J. Blomvall, P.O. Lindberg / European Journal of Operational Research 143 (2002) 452–461
Fig. 4. Convergence of the objective function towards the optimal objective funtion value.
nonlinear separable objective function, 15,777,740 variables and 16,888,850 constraints in 20 hours. The number of scenarios was 1,000,000. Thus we have developed a solver that can solve very large scale multistage stochastic programming problems that have linear or nonlinear separable objective function. The approach presented is quite versatile. In a forthcoming paper, we will give an extension to more general stochastic programming problems. Further, the approach obviously is easily parallelizable. We are in the process to develop a parallel version for a PC-cluster.
We now have to modify the matrices that premultiply the control variables. For the constraint (1b) we make the following modifications:
Appendix A. Initialization
where M is a large positive value. In order to find a starting solution we first initialize ui in a way to make as many as possible of the constraints (1c) feasible. It is only necessary to introduce new control variables, uiP , for those constraints that are not satisfied. This reduces the number of new variables. We want the recently introduced variables, uiP , to become negative, so that they can be eliminated from the optimization problem. To achieve negative uiP the following term is added to the objective function:
We have to initialize the algorithm if we do not have a feasible solution to start from. We want to find a y that satisfies (1c)–(1e) in order to initialize the interior point method. To find a solution that satisfies (1c) we introduce additional control variables, uiP . To simplify the notation uiP has the same dimension as ei . This gives us a new vector of control variables ðu0i ÞT ¼ uTi ; ðuiP ÞT i 2 I:
B0i ¼ ðBi ; 0Þ i 2 I n f1g: In the inequality constraint (1c), we use the variable uiP to find a feasible solution. This leads to the following modifications: Fi0 ¼ ðFi ; I Þ
i 2 I;
where I denotes the identity matrix. Finally (1e) is modified as
0 uli ¼ uli ; 1M i 2 I;
0 uui ¼ uui ; 1M i 2 I;
J. Blomvall, P.O. Lindberg / European Journal of Operational Research 143 (2002) 452–461 1 ðI u0 2l i i
2
gi Þ ;
where 0 0 Ii ¼ 0 I
i2I
gi is the target values for uiP and l is a parameter. References Bertsekas, D.P., 1995. Dynamic Programming and Optimal Control. Athena Scientific, Belmont, MA. Birge, J.R., 1997. Stochastic programming computation and applications. INFORMS Journal on Computing 9 (2), 111– 133. Birge, J., Holmes, D., 1992. Efficient solution of two-stage stochastic linear programs using interior point methods. Computational Optimization and Applications 1 (3), 245–276. Blomvall, J., Lindberg, P.O. Back-testing the performance of an actively managed option portfolio at the Swedish stock
461
market, 1990-1999, Journal of Economic Dynamics and Control, in press. Gondzio, J., Kouwenberg, R., 2001. High performance computing for asset liability management. Operations Research 49, 879–891. Kelly, J.L., 1956. A new interpretation of information rate. Bell System Technical Journal 35, 917–926. Salinger, D.H., Rockafellar, R.T., 1999. Dynamic splitting: An algorithm for deterministic and stochastic multiperiod optimization, working paper, Department of Mathematics, University of Washington, Seattle. Steinbach, M.C., 2000. Hierarchical sparsity in multistage convex stochastic programs, Report 00-15, ZIB, Berlin. Van Slyke, R., Wets, R., 1969. L-shaped linear programs with applications to control and stochastic programming. SIAM Journal on Applied Mathematics 17, 638– 663. Vladimirou, H., Zenios, S.A., 1999. Scalable parallel computations for large-scale stochastic programming. Annals of Operations Research 90, 87–129. von Neumann, J., Morgenstern, O., 1953. Theory of Games and Economic Behavior. Princeton University Press, Princeton.