A Scalable Decomposition Approach for Stochastic Mixed-Integer - PowerPoint PPT Presentation

A Scalable Decomposition Approach for Stochastic Mixed-Integer Programs Kibaek Kim Jointly work with Brian Dandurand (ANL), Cosmin Petra (LLNL), and Victor Zavala (UW-Madison) Mathematics and Computer Science Division Argonne National Laboratory UG Workshop, Berlin, Germany January 15, 2019 1

Argonne: Vital part of DOE National Laboratory System Argonne: Vital part of DOE National Laboratory System 2

Today’s Talk is about … § Stochastic Programming – Integer variables – Parallel algorithms – Software implementation § Dual Decomposition Ax 0 = b 0 T 1 x 0 + W 1 x 1 = b 1 – New parallel algorithms to accelerate solutions . . ... . . . . T N x 0 + W N x N = b N – Asynchronous parallel computation y 0 A + y 1 T 1 + y N T N = π 0 · · · y 1 W 1 = π 1 – New branch-and-bound for global optimality . ... . . + y N W N = π N § Numerical Results – Open-source software package – SIPLIB test instances – Test problems for stochastic unit commitment – Argonne’s high performance computing clusters 3

Parallel Decomposition Methods for Stochastic MIP Kim, Kibaek and Victor M. Zavala. "Algorithmic innovations and software for the dual decomposition method applied to stochastic mixed-integer programs." Mathematical Programming Computation (2017): 1-42. 4

Stochastic Mixed Integer Programming (SMIP) General formulation of SMIP: Recourse function of an integer program: q ( ω ) T y c T x + E [ Q ( x, ω )] Q ( x, ω ) = min min s.t. W ( ω ) y ≥ h ( ω ) − T ( ω ) x, s.t. Ax ≥ b, y ∈ R n 2 − p 2 × Z p 2 x ∈ R n 1 − p 1 × Z p 1 + + + + Make here-and-now decision x Make wait-and-see decision y for given first- • Operational decisions stage decision x and event ω • Logical decisions • Recourse action to event realization • Countable items (e.g. re-scheduling, system restoration) • Time-dependent decisions Observe stochastic event ω • System Failure • Demand and supply • Cost and price • Weather 5

Dual Decomposition for SMIP c T x ( ω 1 ) + q ( ω 1 ) T y ( ω 1 ) + c T x ( ω 1 ) + q ( ω N ) T y ( ω N ) min + . . . − 30 − 35 s.t. Ax ( ω 1 ) ≥ b, c T x + E[Q(x, ω )] − 40 T ( ω 1 ) x ( ω 1 ) + W ( ω 1 ) y ( ω 1 ) ≥ h ( ω 1 ) − 45 . . ... . . − 50 . . − 55 Ax ( ω N ) ≥ b 5 Nonanticipaticity constraints T ( ω N ) x ( ω N ) + W ( ω N ) y ( ω N ) ≥ h ( ω N ) 4 x ( ω 1 ) . . . = x ( ω N ) 3 ∈ R n 1 − p 1 × Z p 1 x ( ω 1 ) , x ( ω N ) . . . 2 + + ∈ R n 2 − p 2 × Z p 2 y ( ω 1 ) , y ( ω N ) . . . 1 5 4 + + x 2 3 2 0 1 0 x 1 Becomes More Challenging! § The recourse function is nonconvex and discontinuous. – Benders Decomposition cannot be used. – Dual Decomposition § Lagrangian relaxation of the nonanticipativity constraints – N c T x j + q T X � � z = min p j j y j x 0 ,x j ,y j N j =1 [ p j ( c T x j + q T j y j ) − λ T H j x j ] X D ( λ ) := min s.t. ( x j , y j ) ∈ G j , j = 1 , . . . , N, x j ,y j j =1 Nonanticipaticity N s.t. ( x j , y j ) ∈ G j , j = 1 , . . . , N X constraints H j x j = 0 ( λ ) j =1 – Seek for the best lower bound by solving N Decomposed X z ≥ z LD := max D j ( λ ) in each scenario j λ j =1 6

Dual-Search Method 1: START Subgradient Method SOLVE the Lagrangian dual D(λ) for given λ k WIDELY USED in non-differentiable § YES D(λ) > z LB ? optimization NO The dual variable is updated as § YES z LB improved for the last m iterations? λ k +1 = λ k − α k X H s x k s NO s ∈ S Subgradient of D( λ ) at λ k REDUCE β k+1 by a half where the step size rule [1] is given by UPDATE z LB z UB � D ( λ k ) α k := β k s k 2 k P s ∈ S H s x k 2 UPDATE λ k+1 and z UB No way of proving optimality § NO min{Duality gap, β k+1 } < ε? Increment k 1. Fisher, Marshall L. "The Lagrangian relaxation method for YES solving integer programming problems." Management science 50.12_supplement (2004): 1861-1871. END 7

Dual-Search Method 2: Cutting Plane Method START Outer approximation to solve § SOLVE the Lagrangian dual D(λ) for given λ 0 X max D s ( λ ) λ ADD linear inequalities to the Master s ∈ S Piecewise linear concave in λ [1] § + λ T ( H s x s ) c T x s + q T � � SOLVE the Master D s ( λ ) = min p s s y s x s ,y s s.t. ( x s , y s ) ∈ conv.hull( G s ) SOLVE the Lagrangian dual D(λ) for given λ k+1 The Master problem is given by § UPDATE z LB if D(λ k+1 ) > z LB X m k = max θ s θ s , λ s ∈ S s ) T ( λ − λ l ) , θ s ≤ D s ( λ l ) + ( H s x l NO s.t. m k- <= D(λ k+1 )? s ∈ S , l = 0 , 1 , . . . , k Increment k YES Finite convergence with Optimality § END 1. Fisher, Marshall L. "The Lagrangian relaxation method for solving integer programming problems." Management science 50.12_supplement (2004): 1861-1871. 8

Dual-Search Method 2: Cutting Plane Method D s ( λ ) D s ( λ ) λ λ 9

New Dual Decomposition Method The cutting-plane solutions of the Master problem can § 15 Subgradient CPM – oscillate significantly; IPM – suffer from degeneracy. 10 || λ k+1 − λ k || 2 Interior-Point Cutting-Plane Method § 5 – Solve the master by using interior-point method to find feasible solutions (vs. extreme optimal points) 0 10 20 30 40 50 60 70 80 90 100 Number of Iterations – Also use the best known upper bound for early termination of IPM iterations Best known upper bound Adding a set of valid inequalities: § – Benders-type feasibility cuts D s ( λ ) – Benders-type optimality cuts – Other cuts can also be added (e.g., cover cuts) – Further cut away the suboptimal solution space λ 10

DSP : Scalable Decomposition Solver Decomposition methods for Structured Programming § Ax 0 = b 0 T 1 x 0 + W 1 x 1 = b 1 . . ... – Exploiting block-angular structures . . . . T N x 0 + W N x N = b N – Dantzig-Wolfe decomposition y 0 A + y 1 T 1 + y N T N = π 0 · · · + (Parallel) Branch-and-Bound y 1 W 1 = π 1 . ... . – Benders decomposition . + y N W N = π N – Dual decompositio Parallel Parallel computing via MPI library § 11

12

DSP reads models from Julia ( 3 ) X min − 1 . 5 x 1 − 4 x 2 + p s Q ( x 1 , x 2 , ξ s 1 , ξ s 2 ) : x 1 , x 2 ∈ { 0 , . . . , 5 } , s =1 where Q ( x 1 , x 2 , ξ s 1 , ξ s 2 ) = min − 16 y 1 + 19 y 2 + 23 y 3 + 28 y 4 y 1 ,y 2 ,y 3 ,y 4 s.t. 2 y 1 + 3 y 2 + 4 y 3 + 5 y 4 ≤ ξ s 1 − x 1 6 y 1 + y 2 + 3 y 3 + 2 y 4 ≤ ξ s 2 − x 2 Only 15 lines of Julia script! y 1 , y 2 , y 3 , y 4 ∈ { 0 , 1 } and ( ξ s 1 , ξ s 2 ) ∈ { (7 , 7) , (11 , 11) , (13 , 13) } with probability 1 / 3. using Dsp, MPI # Load packages 1 MPI.Init() # Initialize MPI 2 m = Model(3) # Create a Model object with three scenarios 3 xi = [[7,7] [11,11] [13,13]] # random parameter 4 @variable (m, 0 <= x[i=1:2] <= 5, Int ) 5 @objective (m, Min, -1.5*x[1]-4*x[2]) 6 for s in 1:3 7 q = Model(m, s, 1/3); 8 @variable (q, y[j=1:4], Bin) 9 @objective (q, Min, -16*y[1]+19*y[2]+23*y[3]+28*y[4]) 10 @constraint (q, 2*y[1]+3*y[2]+4*y[3]+5*y[4]<=xi[1,s]-x[1]) 11 @constraint (q, 6*y[1]+1*y[2]+3*y[3]+2*y[4]<=xi[2,s]-x[2]) 12 end 13 solve(m, solve_type=:Dual, param="myparams.txt") 14 MPI.Finalize() # Finalize MPI 15 13

SIPLIB Test Instances SIPLIB: publically available test library for stochastic integer programming § http://www2.isye.gatech.edu/~sahmed/siplib/ – – DCAP : complete recourse, first-stage mixed -integer and second-stage pure integer • Number of scenarios: 200, 300 and 500 SSLP : complete recourse, first-stage pure integer and second-stage mixed -integer – • Number of scenarios: 5, 10, 15, 50 and 100 Name # Rows # Columns # Integers Name # Rows # Columns # Integers 3006 5412 5406 1501 6505 6255 dcap233 200 sslp 5 25 50 4506 8112 8106 3001 13005 12505 dcap233 300 sslp 5 25 100 7506 13512 13506 301 3465 3390 dcap233 500 sslp 15 45 5 3606 7212 7206 dcap243 200 601 6915 6765 sslp 15 45 10 5406 10812 10806 dcap243 300 901 10365 10135 sslp 15 45 15 9006 18012 18006 dcap243 500 3001 25510 25010 sslp 10 50 50 2406 4812 4806 dcap332 200 6001 51010 50010 sslp 10 50 100 3606 7212 7206 dcap332 300 6006 12012 12006 dcap332 500 2806 6412 6406 dcap342 200 4206 9612 9606 dcap342 300 7006 16012 16006 dcap342 500 14