Global Optimization of ODE constrained network problems Marc Pfetsch Joint work with: Oliver Habeck and Stefan Ulbrich TU Darmstadt Supported by the SFB Transregio 154: Mathematical Modelling, Simulation and Optimization using the Example of Gas Networks March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 1
Outline General Theory 1 Abstract Model and Relaxations Solving Strategies 2 Application – Stationary Gas Transport Model Relaxation of the ODE 3 Implementation – for Stationary Gas Transport The Algorithm Results March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 2
Optimization Problem We consider the ODE constrained optimization problem: C ( x , y 0 , y S , z ) min G ( x , y 0 , y S , z ) ≤ 0, s.t. � � s ∈ [0, S ] ∂ s y ( s ) = f s , x , y ( s ) , ( P ode ) y 0 = y (0), y S = y ( S ), x ∈ X , y 0 ∈ Y 0 , y S ∈ Y S , z ∈ Z , ⊲ Y 0 , Y S ⊂ R n and X ⊂ R k are polytopes and Z ⊂ Z m is bounded. ⊲ C : X × Y 0 × Y S × Z → R and G : X × Y 0 × Y S × Z → R l are continuously differentiable and possibly nonlinear. ⊲ Values of y only need to be known for finitely many points y (0), y ( S ) ∈ R n . ⊲ Examples: Stationary networks (gas, water, (electricity), . . . ) March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 3
Idea of our Approach ⊲ Standard approach is “Discretize then optimize”: discretize ODE (add N variables), use numerical approximation scheme and obtain MINLP . Solve it and possibly refine. ⊲ Our goal: compute globally optimal solution without explicitly discretizing. ⊲ Workhorse: good convex/concave under/over-estimators. ⊲ Will show: Easy to construct in certain cases, e.g., gas optimization. ⊲ Yields globally convergent solution algorithm. March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 4
Equivalent Optimization Problem Equivalent reformulation with the solution mapping F : ( x , y 0 ) �→ y S : C ( x , y 0 , y S , z ) min G ( x , y 0 , y S , z ) ≤ 0, s.t. ( P ) y S − F x , y 0 � � = 0, x ∈ X , y 0 ∈ Y 0 , y S ∈ Y S , z ∈ Z . March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 5
Nonconvex Relaxation of ( P ) Assume that there exist F ℓ : X × Y 0 × N n → R n and F u : X × Y 0 × N n → R n , with F ℓ � x , y 0 , N x , y 0 � ≤ F u � x , y 0 , N � � � ≤ F for all x ∈ X and y 0 ∈ Y 0 . Furthermore, on the polytopes X , Y 0 the functions F ℓ and F u converge uniformly to F for N → ∞ , i = 1, ... , n . This yields the relaxation C ( x , y 0 , y S , z ) min G ( x , y 0 , y S , z ) ≤ 0, s.t. ( P r ( N )) ≤ y S ≤ F u � F ℓ � x , y 0 , N � x , y 0 , N � , x ∈ X , y 0 ∈ Y 0 , y S ∈ Y S , z ∈ Z . March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 6
Convex Relaxation of ( P r (N) ) Assume that there exist F ℓ ≤ F ℓ ⊲ convex underestimators ˇ C ≤ C , ˇ G ≤ G , and ˇ ⊲ and a concave overestimator F u ≤ ˆ F u . Hence, C ( x , y 0 , y S , z ) ˇ min G ( x , y 0 , y S , z ) ≤ 0, ˇ s.t. ( P cv ( N )) ≤ y S ≤ ˆ F ℓ � x , y 0 , N F u � x , y 0 , N ˇ � � , x ∈ X , y 0 ∈ Y 0 , y S ∈ Y S , z ∈ conv( Z ) is a convex relaxation of ( P r ( N )). � we can solve ( P r ( N )) with spatial branch-and-bound. March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 7
( ǫ , δ )-optimal solutions Definition A vector ( x , y 0 , y S , z ) ∈ X × Y 0 × Y S × Z is ( ǫ , δ )-optimal, if ⊲ each constraint is violated by a maximum of δ ⊲ and the objective function satisfies C ( x , y 0 , y S , z ) ≤ C ∗ + ǫ . March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 8
Solving Strategy Goal: Find ( ǫ , δ 1 + δ 2 )-optimal solution of ( P ). ⊲ Start spatial branch-and-bound with some initial N ∈ N n and feasibility tolerance δ 1 . ⊲ Solve the convex relaxation ( P cv ( N )) in a node. ⊲ Check if � F u � x , y 0 , N − F ℓ � x , y 0 , N � � � ∞ ≤ δ 2 holds for the current solution of the relaxation. ⊲ Increase N if necessary, improve the under- and overestimators ˇ F ℓ , ˆ F u , and solve the convex relaxation again. ⊲ Continue like in the “normal” spatial branch-and-bound algorithm. March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 9
Solving Strategy Theorem For a nested sequence of nodes F k = X k × Y 0 k × Y T k × Z k with F k +1 ⊆ F k let ⊲ lim k →∞ diam F k = 0 ⊲ and the estimators become exact in the limit, i.e., the estimators converge to the original functions on F k for k → ∞ . Then the solving strategy produces an ( ǫ , δ 1 + δ 2 ) feasible point of ( P ) or shows that it is infeasible in finite time. March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 10
Outline General Theory 1 Abstract Model and Relaxations Solving Strategies 2 Application – Stationary Gas Transport Model Relaxation of the ODE 3 Implementation – for Stationary Gas Transport The Algorithm Results March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 11
Model ⊲ Gas networks are given by directed graphs G = ( V , A ). ⊲ The nodes are junctions of the network. ⊲ The arcs are (control)valves: linear constraints with binary variables resistors: (non)linear constraints with indicator variables compressors: (here) polyhedral approximation of operating range pipes: ordinary differential equation March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 12
The Differential Equation Stationary isothermal Euler equation with no inclination and | v | c = c | q | Ap ≤ 0.8: λ c 2 q | q | p ( x ) � � ∂ x p ( x ) = − A 2 p ( x ) 2 − c 2 q 2 � =: φ p ( x ), q . � 2 D Variables: ⊲ pressure p ( x ) ⊲ mass flow q Constants: ⊲ friction coefficient λ ⊲ cross sectional area A ⊲ speed of sound c ⊲ length of the pipe L ⊲ diameter D For q ≥ 0 the solution p ( x ) is concave and non increasing. March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 13
Optimization Problem Thus, we want to solve the following optimization problem: min C ( p , q , z ) G ( p , q , z ) ≤ 0, s.t. � � ∀ a ∈ A pipe ⊆ A , ∂ x p a ( x ) = φ a p a ( x ), q a ∀ a = ( u , v ) ∈ A pipe , p u = p a (0), p v = p a ( L a ) q a ≤ q a ≤ q a ∀ a ∈ A , p u ≤ p u ≤ p u ∀ u ∈ V , z ∈ Z ⊂ Z m with spatial branch-and-bound. March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 14
Bounding Methods for Gas Transport Theorem Given a discretization 0 = x N < x N − 1 < · · · < x 0 = L with h i +1 = x i − x i +1 for i = 0, 1, ... , N − 1 . If q ≥ 0 , bounds on p (0) are given by the following methods: ⊲ Upper bounds through the trapezoidal rule: p u p u i +1 + 1 p u = p u i − 1 p u � � � � 0 = p out , 2 h i +1 φ i +1 , q 2 h i +1 φ i , q , ∀ i = 0, 1, ... , N − 1. ⊲ Lower bounds through the explicit midpoint method: p ℓ p ℓ i +1 = p ℓ p ℓ i − 1 p ℓ � � � � 0 = p out , i − h i +1 φ 2 h i +1 φ i , q , q , ∀ i = 0, 1, ... , N − 1. � yields (nonconvex) relaxation of F a . Critical property: nonnegative local truncation error. March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 15
Illustration p in p in F a p in p out p out p out March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 16
Illustration p in p in F a p u ( p out , q , N ) p ℓ ( p out , q , N ) p in p out p out p out March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 16
Illustration p in p in F a R p in p out p out p out March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 16
Illustration p in p in F a conv( R ) p in p out p out p out March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 16
Outline General Theory 1 Abstract Model and Relaxations Solving Strategies 2 Application – Stationary Gas Transport Model Relaxation of the ODE 3 Implementation – for Stationary Gas Transport The Algorithm Results March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 17
Some Informations ⊲ We use the branch-and-bound framework SCIP . ⊲ We implemented a constraint handler for the ODEs. ⊲ Test instances are available at http://gaslib.zib.de/ . March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 18
Basic Algorithm We perform the following steps in the nodes of the branch-and-bound tree: 1. Perform bound propagation (use estimators) 2. Solve LP-relaxation 3. Enforce feasibility of “most infeasible” pipe by one of the following steps: 3.1 Fixing the direction of the flow by branching. 3.2 Adding (parts of) the concave overestimator. 3.3 Branching w.r.t. input pressure, output pressure or massflow. 3.4 Adding a gradient cut (convex underestimator). 4. Resolve new LP-relaxation or choose a new node. March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 19
Illustration p in p in p in p in p in p in p out p out p out p out p out p out March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 20
Example: GasLib-40 ⊲ 40 nodes ⊲ 39 pipes ⊲ 6 compressors ⊲ 12 binary variables for compressors March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 21
Results v ∈V q ± objective: minimize − � v ∈V p v � v p v � a ∈A cs z a solving time (seconds) 877.74 542.48 65.08 processed Nodes 12,444 2790 263 branchings on flow (%) 46.2 81.0 84.2 branchings on pressure (%) 52.9 0.0 8.8 branchings on binary variables (%) 0.9 19.0 7.0 leaves 5900 1412 73 cut offs by propagation 5297 1149 38 bound changes by propagation 32,444 63,317 6051 added overestimators 7362 2385 985 added underestimators 6360 4281 488 δ 1 = ǫ = 10 − 6 , δ 2 = 10 − 4 March 7, 2018 | SCIP Workshop 2018 | M. Pfetsch | 22
Recommend
More recommend
