Convex Relaxations in Mixed-Integer Optimization Methods and - - PowerPoint PPT Presentation

Convex Relaxations in Mixed-Integer Optimization – Methods and Control Applications

Robin Vujanic

Institut f¨ ur Automatik (IfA) Department of Electrical Engineering Swiss Federal Institute of Technology (ETHZ)

August 20, 2014

1 / 19

Outline

1 Introduction 2 Optimization of Large Scale Systems 3 Example: Electric Vehicles Charging Coordination

2 / 19

What are Mixed–Integer Optimization Problems / Why do we look into them

• many practical and industrial systems entail continuous quantities

◮ physical measurements of voltages ◮ concentrations ◮ and positions in space

as well as discrete components

◮ on/off decisions ◮ switches ◮ and logic reasoning (if, or, . . . )

• when the associated control/operation tasks are tackled with
• ptimization, Mixed-Integer Optimization Problems (MIPs) arise
• but the additional flexibility has a price ...

3 / 19

Computational Issues – a Practical Perspective

• experience with instances from supply chain problem
• modest size, yet memory blew up

# cust. # prod. CPLEX time (sec) 100 25 Min: 10.3 Avg: 63.5 Max: 202.5 300 75

• ut of memory
• computational requirements depend not only on structure and size of

the problem, but also the data

◮ bad in a control context

The thesis focuses on particular model structures that are of practical

• interest. For these, we derived computationally attractive

approximation schemes, equipped with guarantees.

4 / 19

Thesis Content

Part I – Optimization of Large Scale Systems.

• Lagrangian Duality
• Applications:

◮ Power Systems ◮ Supply Chain Optimization

Part II – Robust Optimization of Uncertain Systems.

• Linear Programming Relaxations
• Applications:

◮ Scheduling under Uncertainty ◮ PWM Systems 5 / 19

• [EEM-13]
• [MED-14]
• [CDC-14]
• [Submitted, Math. Prog.-13]
• [ACC-12]
• [ECC-13]
• [CDC-13]

Thesis Content

Part I – Optimization of Large Scale Systems.

• Lagrangian Duality
• Applications:

◮ Power Systems ◮ Supply Chain Optimization

Part II – Robust Optimization of Uncertain Systems.

• Linear Programming Relaxations
• Applications:

◮ Scheduling under Uncertainty ◮ PWM Systems 5 / 19

• [EEM-13]
• [MED-14]
• [CDC-14]
• [Submitted, Math. Prog.-13]
• [ACC-12]
• [ECC-13]
• [CDC-13]

Outline

1 Introduction 2 Optimization of Large Scale Systems 3 Example: Electric Vehicles Charging Coordination

6 / 19

Large Scale Optimization – Problem Structure Considered

• we consider the problem

P :    min

• i cixi

s.t.

• i Hixi ≤ b

xi ∈ Xi i = 1, . . . , I where Xi = {xi ∈ Rri × Zzi | Aixi ≤ di}

• large collection of subsystems

◮ subsystem model is Xi (mixed–integer) ◮ coupled by shared resources → coupling constraints

i Hixi ≤ b

• # of subsystems I ≫ # of coupling constraints m

7 / 19 . . .

Substation

. . . . . . . . . . . . . . .

aggregated flow Rest of the Network critical branch flow

Problem’s Decomposition

Obtain decomposition using duality: minx

• i∈I cixi

s.t.

• i∈I Hixi ≤ b

xi ∈ Xi ⇒ minx

• i∈I cixi + λ′(

i∈I Hixi − b)

s.t. xi ∈ Xi ⇒

• i∈I

min

xi∈Xi

• cixi + λ′(Hixi)
• − λ′b
• .

=d(λ)

• Lagrangian dual (or outer) problem:

D : max d(λ) s.t. λ ≥ 0

8 / 19

Solutions to the inner problem

d(λ) =

• i∈I

min

xi∈Xi

• cixi + λ′(Hixi)
• − λ′b

Consider xi(λ⋆) ∈ arg min

xi∈Xi

• cixi + λ⋆(Hixi)
• as candidate solution to P

Properties of xi(λ⋆):

• satisfy Xi constraints
• obtained “for free” as by-product of methods that solve D
• distributed computations
• generally infeasible in the MIP case !

◮ violate coupling constraints 9 / 19

Primal Recovery Scheme

• we show that in x(λ⋆) only m subsystems may be “problematic”

◮ technique based on Shapley–Folkman–Starr theorem [CDC ’14] ◮ or simplex tableaux argument [MED ’13] ◮ these arguments are used to show bounded duality gap [Ekeland ’76,

Bertsekas ’83]

• so we propose to consider instead

P :    minx

• i∈I cixi

s.t.

• i∈I Hixi ≤ b − ρ

xi ∈ Xi ∀i ∈ I, where ρ = m · max

i∈I

• max

xi∈XiHixi − min xi∈XiHixi

• Theorem

Then x(¯ λ⋆) is feasible for P. [under some uniqueness assumptions]

10 / 19

Performance of the Recovered Solutions

• under some technical assumption, ...

Theorem

The recovered solution x(¯ λ⋆) is feasible and satisfies JP(x(¯ λ⋆)) − J⋆

P ≤ (m + ||ρ||∞/ζ) ·

• max

xi∈Xi cixi − min xi∈Xi cixi

• if J⋆

P grows linearly with |I|, and Xi uniformly bounded

J(x(¯ λ⋆)) − J⋆

P

J⋆

P

→ 0 as |I| → ∞

11 / 19

Improving Approximations – Conservatism Reduction

ρ scales with m but not with I – want to keep it as small as possible

• when couplings are determined by certain network topologies

A B C D E F 1 2 3 4 5 6 7 8 9 10 11 12 13 15 16 17 18 19 14 20 21 22 23 25 26 27 24 28 Subsystems A B C D E F 7 13 14 19 23 20 24 28 − − − − k-th row [Hi]i∈Ik 1 6 − A1 H1 H2 HI A2 AI . . . . . . . . . . . . . . . . . .

• can safely use rank([Hi]i∈Ik) instead of m
• generally possible to use rank(H) instead of m

12 / 19

Outline

1 Introduction 2 Optimization of Large Scale Systems 3 Example: Electric Vehicles Charging Coordination

13 / 19

Electric Vehicle (EV) Charging Coordination [CDC ’14]

• expected increase of EV presence
• substantial additional stress on

network & equipment ⇒ need charging coordination

• network administrator (DSO)

can’t tackle each single unit ⇒ EV aggregator

14 / 19

Aggregator’s Role

• aggregator’s control task is to

assign to each EV the time slots when charging can occur

• compatibly with...
• local requirements

◮ required final State of Charge ◮ fixed charge rates ◮ battery capacity limits

• global objectives

◮ network congestion avoidance

(limits set by DSO)

◮ “valley fill”, cost min., . . . 15 / 19

from [Lopes ’11] 0:00 1:30 3:00 4:30 6:00 7:30 10 20 State of Charge (kWh) SoC desired final SoC 0:00 1:30 3:00 4:30 6:00 7:30 0.5 1 Charge Control Time of the Day

Computational Experiments

• cast as large optimization problem
• solve using proposed method: duality+contraction

◮ support extensions (e.g., vehicle–to–grid “V2G”)

• population up to 10’000 EVs
• computation times ≤ 10 sec (charge only)

◮ greedy subproblem structure 16 / 19

Solutions – Charge and V2G

0:00 1:30 3:00 4:30 6:00 7:30 2 4 6 8 10 Time of the Day Power Profiles (MW) reference iteration #80 final iteration

(a) reference tracking

14:00 17:00 20:00 23:00 2:00 5:00 8:00 11:00 10 20 30 40 50 60 70 Time of the Day Total Power Flow (MW) base load base + EVs load

(b) resulting “valley fill”

0:00 1:30 3:00 4:30 6:00 7:30 0.5 1 1.5 2 2.5 3 Time of the Day Branch Power Flow (MW) line capacity contracted line cap. flow iteration #80 final iteration

(c) network limits

0:00 1:30 3:00 4:30 6:00 7:30 10 20 State of Charge (kWh) SoC desired final SoC 0:00 1:30 3:00 4:30 6:00 7:30 −1 1 Charge Control Time of the Day

(d) local requirements

17 / 19

Other Examples or Applications

• supply chains optimization – partial shipments [MED ’14]
• power systems operation

◮ control of TCLs ◮ large fleet of generators

• portfolio optimization for small investors
• . . .

18 / 19

Questions?

• Prof. Manfred Morari
• Dr. Paul Goulart
• Dr. S´

ebastien Mari´ ethoz

• Dr. Peyman Mohajerin

Esfahani

• Dr. Joe Warrington
• Dr. Stefan Richter
• Dr. Sumit Mitra
• Gregory Ledva
• Isik Ilber Sirmatel
• Marius Schmitt
• Prof. Debasish Chatterjee
• Prof. Federico Ramponi
• Dr. Peter Hokayem
• Dr. Apostolos Fertis
• Dr. Utz-Uwe Haus
• Dr. Alexander Fuchs
• Dr. Martin Herceg
• Robert Nguyen

19 / 19

References

• [Ekeland ’76]J. P. Aubin and I. Ekeland, Estimates of the duality gap

in nonconvex optimization, Mathematics of Operations Research 1 (1976), no. 3, 225-–245.

• [Bertsekas ’83] Dimitri P. Bertsekas, G. Lauer, N. Sandell, and T.

Posbergh, Optimal short-term scheduling of large-scale power systems, IEEE Transactions on Automatic Control 28 (1983), no. 1, 1– 11.

• [Dawande ’06] Milind Dawande, Srinagesh Gavirneni, and Sridhar

Tayur, Effective heuristics for multiproduct partial shipment models, Operations Research 54 (2006), no. 2, 337–352 (en).

• [Vujanic ’13] R. Vujanic, P. Mohajerin Esfahani, P. Goulart, S.

Mariethoz and M. Morari, Vanishing Duality Gap in Large Scale Mixed-Integer Optimization: a Solution Method with Power System Applications, submitted to Mathematical Programming (2013).

• [Lopes ’11] J. Lopes, F. Soares, and P. Almeida, Integration of

Electric Vehicles in the Electric Power System, Proceedings of the IEEE, 2011, p.168-183.

19 / 19

BACKUP SLIDES

19 / 19

Convex Relaxations of MIPs

• our methods are based on convex relaxations
• given pure binary problem

P :        minx c⊺x s.t. a1x ≤ b1 a2x ≤ b2 x ∈ {0, 1}n

• if you had to solve it, which relaxation would you pick?

P1 :        minx c⊺x s.t. a1x ≤ b1 a2x ≤ b2 conv(x ∈ {0, 1}n) P2 :        minx c⊺x s.t. a1x ≤ b1 conv a2x ≤ b2 x ∈ {0, 1}n

• r

P3 :        minx c⊺x s.t. conv   a1x ≤ b1 a2x ≤ b2 x ∈ {0, 1}n  

19 / 19

Strength of the Relaxations

X1 X2 min .

• i∈I cixi

s.t.

• i∈I Hixi ≤ b

Aixi ≤ di xi ∈ {0, 1}ni X1 X2 min .

• i∈I cixi

s.t.

• i∈I Hixi ≤ b

Aixi ≤ di conv

• xi ∈ {0, 1}ni

X1 X2 min .

• i∈I cixi

s.t.

• i∈I Hixi ≤ b

conv

  Aixi ≤ di

xi ∈ {0, 1}ni

 

X1 X2 min .

• i∈I cixi

s.t. conv

     

• i∈I Hixi ≤ b

Aixi ≤ di xi ∈ {0, 1}ni

     

19 / 19

Issues affecting inner solutions

x1 x2 1 1 f conv(Xi)

min −x1 s.t. x1 − x2 ≤ 0.5 x1 + x2 ≤ 1.5 x1, x2 ∈ {0, 1}

19 / 19

Issues affecting inner solutions

x1 x2 1 1 f conv(Xi)

min −x1 s.t. x1 − x2 ≤ 0.5 x1 + x2 ≤ 1.5 x1, x2 ∈ {0, 1}

x⋆

19 / 19

Issues affecting inner solutions

x1 x2 1 1 f conv(Xi)

min −x1 s.t. x1 − x2 ≤ 0.5 x1 + x2 ≤ 1.5 x1, x2 ∈ {0, 1}

x⋆ ⇒ max

λ≥0

min

x1,x2∈{0,1}

−x1 + λ1(x1 − x2 − 0.5) + λ2(x1 + x2 − 1.5)

19 / 19

Issues affecting inner solutions

⇒ d⋆ = −1 λ⋆ = (0.25, 0.25) x1 x2 1 1 f conv(Xi)

min −x1 s.t. x1 − x2 ≤ 0.5 x1 + x2 ≤ 1.5 x1, x2 ∈ {0, 1}

x⋆ ⇒ max

λ≥0

min

x1,x2∈{0,1}

−x1 + λ1(x1 − x2 − 0.5) + λ2(x1 + x2 − 1.5)

19 / 19

Issues affecting inner solutions

⇒ d⋆ = −1 λ⋆ = (0.25, 0.25) x1 x2 1 1 f conv(Xi)

min −x1 s.t. x1 − x2 ≤ 0.5 x1 + x2 ≤ 1.5 x1, x2 ∈ {0, 1}

x⋆ ⇒ max

λ≥0

min

x1,x2∈{0,1}

−x1 + λ1(x1 − x2 − 0.5) + λ2(x1 + x2 − 1.5) x(λ⋆)

19 / 19

Issues affecting inner solutions

⇒ d⋆ = −1 λ⋆ = (0.25, 0.25) x1 x2 1 1 f conv(Xi)

min −x1 s.t. x1 − x2 ≤ 0.5 x1 + x2 ≤ 1.5 x1, x2 ∈ {0, 1}

x⋆ ⇒ max

λ≥0

min

x1,x2∈{0,1}

−x1 + λ1(x1 − x2 − 0.5) + λ2(x1 + x2 − 1.5) x(λ⋆) x⋆

LP

19 / 19

Illustration of the Problems with the L-Relaxation

x1 x2 3x1 − x2 x1 − x2 ≥ −1 − x1 + 2x2 ≤ 5 3x1 + x2 ≥ 3 6x1 + x2 ≤ 15 x1, x2 ≥ 0 x1, x2 ∈ Z minx s.t. 1 2 3 1 2 3

19 / 19

Illustration of the Problems with the L-Relaxation

x1 x2 3x1 − x2 x1 − x2 ≥ −1 − x1 + 2x2 ≤ 5 3x1 + x2 ≥ 3 6x1 + x2 ≤ 15 x1, x2 ≥ 0 x1, x2 ∈ Z minx s.t. 1 2 3 1 2 3 x∗

19 / 19

Illustration of the Problems with the L-Relaxation

x1 x2 3x1 − x2 x1 − x2 ≥ −1 − x1 + 2x2 ≤ 5 3x1 + x2 ≥ 3 6x1 + x2 ≤ 15 x1, x2 ≥ 0 x1, x2 ∈ Z minx s.t. 1 2 3 1 2 3 x∗

19 / 19

Illustration of the Problems with the L-Relaxation

x1 x2 3x1 − x2 x1 − x2 ≥ −1 − x1 + 2x2 ≤ 5 3x1 + x2 ≥ 3 6x1 + x2 ≤ 15 x1, x2 ≥ 0 x1, x2 ∈ Z minx s.t. 1 2 3 1 2 3 x∗

19 / 19

Illustration of the Problems with the L-Relaxation

x1 x2 3x1 − x2 x1 − x2 ≥ −1 − x1 + 2x2 ≤ 5 3x1 + x2 ≥ 3 6x1 + x2 ≤ 15 x1, x2 ≥ 0 x1, x2 ∈ Z minx s.t. 1 2 3 1 2 3 x∗ d(λ) = min

x∈X

• 3x1 − x2 + λ(−1 − x1 + x2)
• 19 / 19

Illustration of the Problems with the L-Relaxation

x1 x2 3x1 − x2 x1 − x2 ≥ −1 − x1 + 2x2 ≤ 5 3x1 + x2 ≥ 3 6x1 + x2 ≤ 15 x1, x2 ≥ 0 x1, x2 ∈ Z minx s.t. 1 2 3 1 2 3 x∗ d(λ) = min

x∈X

• 3x1 − x2 + λ(−1 − x1 + x2)
• ⇒

max

λ≥0 d(λ) 19 / 19

Illustration of the Problems with the L-Relaxation

x1 x2 3x1 − x2 x1 − x2 ≥ −1 − x1 + 2x2 ≤ 5 3x1 + x2 ≥ 3 6x1 + x2 ≤ 15 x1, x2 ≥ 0 x1, x2 ∈ Z minx s.t. 1 2 3 1 2 3 x∗ d(λ) = min

x∈X

• 3x1 − x2 + λ(−1 − x1 + x2)
• ⇒

max

λ≥0 d(λ)

⇒ λ∗ = 5/3

19 / 19

Illustration of the Problems with the L-Relaxation

x1 x2 3x1 − x2 x1 − x2 ≥ −1 − x1 + 2x2 ≤ 5 3x1 + x2 ≥ 3 6x1 + x2 ≤ 15 x1, x2 ≥ 0 x1, x2 ∈ Z minx s.t. 1 2 3 1 2 3 x∗ d(λ) = min

x∈X

• 3x1 − x2 + λ(−1 − x1 + x2)
• ⇒

max

λ≥0 d(λ)

⇒ λ∗ = 5/3 ⇒ X(λ∗) = 1

• ,

2

• 19 / 19

Illustration of the Problems with the L-Relaxation

x∗ x1 x2 3x1 − x2 x1 − x2 ≥ −1 − x1 + 2x2 ≤ 5 3x1 + x2 ≥ 3 6x1 + x2 ≤ 15 x1, x2 ≥ 0 x1, x2 ∈ Z minx s.t. 1 2 3 1 2 3 d(λ) = min

x∈X

• 3x1 − x2 + λ(−1 − x1 + x2)
• ⇒

max

λ≥0 d(λ)

⇒ λ∗ = 5/3 ⇒ X(λ∗) = 1

• ,

2

• conv(Xi)

19 / 19

Illustration of the Problems with the L-Relaxation

x∗ x1 x2 3x1 − x2 x1 − x2 ≥ −1 − x1 + 2x2 ≤ 5 3x1 + x2 ≥ 3 6x1 + x2 ≤ 15 x1, x2 ≥ 0 x1, x2 ∈ Z minx s.t. 1 2 3 1 2 3 d(λ) = min

x∈X

• 3x1 − x2 + λ(−1 − x1 + x2)
• ⇒

max

λ≥0 d(λ)

⇒ λ∗ = 5/3 ⇒ X(λ∗) = 1

• ,

2

• conv(Xi)

x(λ∗)

19 / 19

Illustration of the Problems with the L-Relaxation

x∗ x1 x2 3x1 − x2 x1 − x2 ≥ −1 − x1 + 2x2 ≤ 5 3x1 + x2 ≥ 3 6x1 + x2 ≤ 15 x1, x2 ≥ 0 x1, x2 ∈ Z minx s.t. 1 2 3 1 2 3 d(λ) = min

x∈X

• 3x1 − x2 + λ(−1 − x1 + x2)
• ⇒

max

λ≥0 d(λ)

⇒ λ∗ = 5/3 ⇒ X(λ∗) = 1

• ,

2

• conv(Xi)

x(λ∗) ”p∗”

19 / 19

Main Idea

The following convexification of P is strongly related to D (PLP) :    minx

• i∈I cixi

s.t.

• i∈I Hixi ≤ b

xi ∈ conv(Xi) ∀i ∈ I Known that:

• J⋆

PLP = J⋆ D[Geoffrion ’74]

We show that:

1 an optimal solution x⋆ LP to PLP has nice properties

◮ feasible w.r.t. the complicating constraints ◮ returns integral solutions for |I| − m subproblems

2 x⋆ LP and x(λ⋆) differ in at most m subproblems 3 thus we need limited compensation to make x(λ⋆) feasible

19 / 19

Main Idea

The following convexification of P is strongly related to D (PLP) :    minx

• i∈I cixi

s.t.

• i∈I Hixi ≤ b

xi ∈ conv(Xi) ∀i ∈ I Known that:

• J⋆

PLP = J⋆ D (surprising!)

We show that:

1 an optimal solution x⋆ LP to PLP has nice properties

◮ feasible w.r.t. the complicating constraints ◮ returns integral solutions for |I| − m subproblems

2 x⋆ LP and x(λ⋆) differ in at most m subproblems 3 thus we need limited compensation to make x(λ⋆) feasible

19 / 19

The Shapley-Folkman Theorem

Theorem (Shapley-Folkman)

Let Si, i = 1, . . . , |I| be nonempty subsets of Rm, with |I| > m, and let S = S1 + · · · + S|I|. Then every vector s ∈ conv(S) can be represented as s = s1 + · · · + s|I|, where si ∈ conv(Si) for all i = 1, . . . , |I|, and si / ∈ Si for at most m indices i.

• relatively known in economics
• Lloyd Shapley won the Nobel Prize in Economic Sciences (in ’12)
• in our case

Si = HiXi S = ⊕i∈ISi =

• i∈I Hixi | xi ∈ Xi
• 19 / 19

Visualization of the Theorem

minx

• i∈I cixi

s.t.

• i∈I Hixi ≤ b

xi ∈ Xi ∀i ∈ I

0.5 1 0.5 1 X1 0.5 1 X2 0.5 1 0 0.5 1 2 X3 0.5 1 0.5 1 X4

(a) subsystems Xi

−2 2 −3 −2 −1 1 2 3

1 2

○

1 2

19 / 19

Visualization of the Theorem

minx

• i∈I cixi

s.t.

• i∈I Hixi ≤ b

xi ∈ Xi ∀i ∈ I

0.5 1 0.5 1 X1 0.5 1 X2 0.5 1 0 0.5 1 2 X3 0.5 1 0.5 1 X4

(a) subsystems Xi

−2 2 1 2 3 S1 = H1X1 0.5 1 −3 −2 −1 S2 = H2X2 −1 1 −1 1 2 S3 = H3X3 1 2 0.5 1 1.5 S4 = H4X4

(b) budget consumption HiXi

−2 2 −3 −2 −1 1 2 3

1 2

○

1 2

19 / 19

Visualization of the Theorem

minx

• i∈I cixi

s.t.

• i∈I Hixi ≤ b

xi ∈ Xi ∀i ∈ I

0.5 1 0.5 1 X1 0.5 1 X2 0.5 1 0 0.5 1 2 X3 0.5 1 0.5 1 X4

(a) subsystems Xi

−2 2 1 2 3 S1 = H1X1 0.5 1 −3 −2 −1 S2 = H2X2 −1 1 −1 1 2 S3 = H3X3 1 2 0.5 1 1.5 S4 = H4X4

(b) budget consumption HiXi

−2 2 4 6 −3 −2 −1 1 2 3 4 5 6 H1(x) H2(x) S =

iHiXi

○ + b1 b2

(c) aggregated budget con- sumption

19 / 19

Visualization of the Theorem

minx

• i∈I cixi

s.t.

• i∈I Hixi ≤ b

xi ∈ conv(Xi) ∀i ∈ I

0.5 1 0.5 1 X1 0.5 1 X2 0.5 1 0 0.5 1 2 X3 0.5 1 0.5 1 X4

(a) subsystems Xi

−2 2 1 2 3 S1 = H1X1 0.5 1 −3 −2 −1 S2 = H2X2 −1 1 −1 1 2 S3 = H3X3 1 2 0.5 1 1.5 S4 = H4X4

(b) budget consumption HiXi

−2 2 4 6 −3 −2 −1 1 2 3 4 5 6 H 1(x) H 2(x) S and conv(S) b1 b2

(c) aggregated budget con- sumption

19 / 19

Visualization of the Theorem

minx

• i∈I cixi

s.t.

• i∈I Hixi ≤ b

xi ∈ conv(Xi) ∀i ∈ I

0.5 1 0.5 1 X1 0.5 1 X2 0.5 1 0 0.5 1 2 X3 0.5 1 0.5 1 X4

(a) subsystems Xi

−2 2 1 2 3 S1 = H1X1 0.5 1 −3 −2 −1 S2 = H2X2 −1 1 −1 1 2 S3 = H3X3 1 2 0.5 1 1.5 S4 = H4X4

(b) budget consumption HiXi

−2 2 4 6 −3 −2 −1 1 2 3 4 5 6 H 1(x) H 2(x) S and conv(S) b1 b2

(c) aggregated budget con- sumption

19 / 19

Result #1: Bounded Duality Gap

• Duality gap for programs structured as P is bounded:

Theorem (bound on duality gap (Ekeland ’76, Bertsekas ’83))

Assume that the sets Xi are non-empty and compact, and that for any given xi ∈ conv(Xi), there exists an ˜ xi ∈ Xi such that Hi ˜ xi ≤ Hixi. Then J⋆

P − J⋆ D ≤ m · max i∈I

• max

xi∈Xicixi − min xi∈Xicixi

• non-convexities do not compound each other indefinitely in P, i.e. the

distance between P and PLP does not grow

• Thus, if J⋆

P increases linearly with |I|,

J⋆

P − J⋆ D

J⋆

P

→ 0 as |I| → ∞.

19 / 19

Example #2: Partial Shipments

• distributor has to supply I customers
• often, available inventory < total demand

◮ demand uncertainty and high storage costs ◮ restrictions on manufacturing capacity

• thus

◮ EITHER fully supply a smaller set of customers ◮ OR partially supply a larger set

• problem is to decide which customer gets what product in the

presence of shipping restrictions

19 / 19

Example #2: Partial Shipments

• distributor has to supply I customers
• often, available inventory < total demand

◮ demand uncertainty and high storage costs ◮ restrictions on manufacturing capacity

• thus

◮ EITHER fully supply a smaller set of customers ◮ OR partially supply a larger set

• problem is to decide which customer gets what product in the

presence of shipping restrictions

19 / 19

Optimization Problem Model

i = customer, j = product wi = ship to i (yes/no), Sj

i = fraction shipped

                   max

• i

reward for wi +

i,j

reward for Sj

i

s.t.

• i

Sj

i ≤ inventory of prod. j

     if wi = 1 ⇒

• j Sj

i ≥ min. shipment

and Sj

i ≤ demand of j

if wi = 0 ⇒ Sj

i = 0

    

• compare with generic model P:

P :    minx

• i∈I cixi

s.t.

• i∈I Hixi ≤ b

xi ∈ Xi ∀i ∈ I

19 / 19

Results using Proposed Method

• partial shipment problem is NP-Hard
• greedy strategies perform poorly
• tailored MINTO setup also not good (see [Dawande ’06])

◮ at N = 300, M = 75

avg opt. gap = 6.2% and computation time 6h

• generate instances of industrial size
• solve with adapted method

Proposed Method CPLEX cust. prod. gap (%)† time (s)† time (s) 100 25 0.38 18.5 Min: 10.3 Avg: 63.5 Max: 202.5 300 75 0.03 62.5 – 500 50 0.06 101.0 – 600 50 0.06 116.6 – 1000 100 0.03 245.2 – † average over 10 instances

19 / 19

x⋆

LP can be computed

1 solve D using a subgradient method

λ[1] = λ(k+1) = P+

• λ(k) + s(k) · γ(k)

(1)

◮ s(k) steplength, pick s(k) = α/k. ◮ γ(k) a subgradient, pick γ(k) =

i∈I Hixi(λ(k)) − b.

2 while iterating (1), calculate average

¯ x(k) = 1 k

k

• j=1

x(λ(j))

Theorem (primal + dual convergence)

• λ(k) → λ⋆ ∈ Λ⋆,
• if x⋆

LP is the unique solution to PLP, ¯

x(k) → x⋆

LP.

19 / 19

Computing approximate solutions to P

Assume Hixi ≥ 0, and 0 ∈ Xi (relaxed in [Vujanic ’13]).

• choices of xi always “consume” budget

Procedure to compute approx. ˆ x⋆ to P (distributed)

1 compute x⋆ LP using averaging 2 find I1 ⊆ I such that (x⋆ LP)i ∈ vert(Xi)

◮ |I1| ≥ |I| − m − 1

3 for i ∈ I1, set ˆ

x⋆

i = (x⋆ LP)i 4 for the remaining i ∈ I \ I1,

◮ EITHER ˆ

x⋆

i = 0

◮ OR set b → b −

i∈I1 Hi ˆ

x⋆

i , and solve smaller dimensional problem

Then, ˆ x⋆ is feasible for P, and satisfies JP(x⋆

ˆ I ) − J⋆ P ≤ (m + 1) max i∈I max xi∈Xi c⊺ i xi.

19 / 19

Control Strategy

• Optimal Control + Lagrangian Duality
• Optimization Problem Model

minimize

e,u,v

• i∈I

Pi (C u · ui − C v · vi) subject to Pmin ≤

i∈I

Pi(ui − vi) ≤ Pmax (ei, ui, vi) ∈ Xi with local (private) model Xi =          ei ui vi   ∈ RN × Z2N

• ei = E init

i

+ B

• E in

i ui − E out i

vi

• E min

i

≤ ei ≤ E max

i

ei[N] ≥ E ref

i

0 ≤ ui + vi ≤ 1       

• ui[k], vi[k] ∈ {0, 1} : charge/discharge decision at step k for EV i ∈ I
• Pi: charge rate
• C u[k], C v[k]: prices for charging/discharging

19 / 19

Solution Method

1 set

ρcharge = N · maxi∈I Pi ρV2G = 2N · maxi∈I Pi

• ⇒

¯ Pmin = Pmin + ρ ¯ Pmax = Pmax − ρ

2 add δu i [k], δv i [k] perturbations to cost function 3 dual (outer) problem

◮ solved using a subgradient method ◮ constant stepsize, decreased every 20–30 iterations

4 inner problem

◮ charge only: optimal local solution is greedy (easy to show) ◮ charge+V2G: deploy DP (local problems are 1D) 19 / 19

Results (5000 EVs) (1/2)

2 4 6 8 −15 −10 −5 5 10 15 time (h)

• i∈I pi: power flow (MW)

Iteration # 20

19 / 19

Results (5000 EVs) (1/2)

2 4 6 8 −15 −10 −5 5 10 15 time (h)

• i∈I pi: power flow (MW)

Iteration # 40

19 / 19

Results (5000 EVs) (1/2)

2 4 6 8 −15 −10 −5 5 10 15 time (h)

• i∈I pi: power flow (MW)

Iteration # 70

19 / 19

Results (5000 EVs) (1/2)

2 4 6 8 −15 −10 −5 5 10 15 time (h)

• i∈I pi: power flow (MW)

Iteration # 100

19 / 19

Results (5000 EVs) (1/2)

2 4 6 8 −15 −10 −5 5 10 15 time (h)

• i∈I pi: power flow (MW)

Iteration # 150 (last)

19 / 19

Results (5000 EVs) (2/2)

50 100 150 1.5 2 2.5 3 3.5 4 4.5 5 iteration # dual objective (k euro) 50 100 150 1 2 3 4 5 iteration # dual objective (k euro) 50 100 150 20 40 60 80 100 iteration #

• (

i∈I Hixi − b)+

• 1 (MW)

lower bound infeas. upper bound infeas.

(c) with disturbances

50 100 150 20 40 60 80 100 iteration #

• (

i∈I Hixi − b)+

• 1 (MW)

lower bound infeas. upper bound infeas.

(d) without disturbances

19 / 19

Results: Charge Only

Proposed Method CPLEX

• Opt. Gap (%)

Solve time (sec) Solve time (sec) # of EVs Min Avg Max Min Avg Max Min Avg Max 200 3.24 3.32 3.41 * * * 1.97 2.16 3.74 350 2.21 2.44 2.58 * * * 1.13 1.79 2.30 500 1.40 1.46 1.54 * * * 1.02 1.24 1.52 700 1.01 1.05 1.10 0.31 0.31 0.31 1.27 1.29 1.31 1000 0.68 0.72 0.76 0.44 0.44 0.44 1.68 1.70 1.73 1500 0.46 0.47 0.49 0.67 0.70 0.70 2.39 2.42 2.45 2000 0.33 0.35 0.36 0.88 0.88 0.89 3.22 3.30 3.41 5000 0.13 0.14 0.14 2.17 2.21 2.23 8.00 8.18 8.43 7000 0.05 0.05 0.06 3.14 3.15 3.16 11.40 11.74 13.25 10000 0.03 0.03 0.04 4.45 4.51 4.52 17.41 17.77 18.43 (*) ≤ 0.3 sec (imprecise measurements).

• all solutions recovered are feasible
• fast computation times due to greedy subproblem structure

19 / 19

Results: Charge + V2G

Proposed Method CPLEX

• Opt. Gap (%)

Solve time (min) Solve time (min) # of EVs Min Avg Max Min Avg Max Min Avg Max 200 8.82 10.51 12.37 1.05 1.06 1.08 0.06 0.07 0.07 350 2.93 3.24 3.51 1.48 1.49 1.52 1.56 6.89 15.81 500 2.05 2.15 2.24 1.85 1.93 2.62 15.21(*) 65.10(*) 262.81(*) 700 1.48 1.54 1.61 2.43 2.44 2.48 – – – 1000 1.01 1.05 1.10 3.24 3.26 3.28 – – – 1500 0.65 0.68 0.72 4.72 4.74 4.81 – – – 2000 0.45 0.50 0.53 6.19 6.21 6.23 – – – 5000 0.12 0.15 0.20 14.88 14.90 14.95 – – – 7000 0.09 0.10 0.12 20.59 20.77 21.94 – – – 10000 0.06 0.07 0.07 29.34 29.39 29.58 – – – (*) failed to solve two instances (out of memory) (–) out of memory before attaining the desired optimality gap

• 500 EVs: CPLEX times vary from 15 min to 4.5h (6+ when unsolved)
• 10k EVs: with 500k binaries and 250k continuous vars, these are among the

largest MIPs publicly available

• CPLEX does not even finish preprocessing on those

19 / 19

Solutions Performance

2000 4000 6000 8000 10000 2 4 6 8 10 12 14 # of EVs

• ptimality gap (%) (avg, min, max)
• “1/|I|” rate of decrease of optimality gap verified in both cases

19 / 19

Intuitive explanation of proof that 1) x⋆

LP is structured and

2) that x⋆

LP and xi(λ) are connected

• PLP is equivalent to

   min

• i cixi

s.t.

• i Hixi ≤ b

xi ∈ conv(Xi) ⇔          minp

• i∈I
• j∈Ji pij(cixj

i )

s.t.

• i∈I
• j∈Ji pij(Hixj

i ) ≤ b

• j∈Ji pj

i = 1

i ∈ I pj

i ≥ 0

i ∈ I, j ∈ Ji

• at most |I| + m entries of p are > 0

◮ p is dimension K (huge), hence requires K equalities to be determined ◮ in PLP, only |I| + m equalities, so K − (|I| + m) equalities have to be

picked form pj

i ≥ 0

◮ hence at most |I| + m entries pj

i can be > 0 (“at most” comes from

slack for coupling constr.)

• for at least |I| − m subproblems (p⋆)ˆ

i i = 1

• inner problem can return x j

i only if its “probability” is (p⋆)j i > 0

• we show this using strict complementarity, which requires

uniqueness

19 / 19

Robust Schedules

Is robustification of the schedule equivalent to, e.g., increasing task length?

19 / 19

Robust Schedules

Is robustification of the schedule equivalent to, e.g., increasing task length?

encode possible delay

• f 1 unit by assuming

task A is length 2

19 / 19

Robust Schedules

Is robustification of the schedule equivalent to, e.g., increasing task length?

encode possible delay

• f 1 unit by assuming

task A is length 2

19 / 19

Robust Schedules

Is robustification of the schedule equivalent to, e.g., increasing task length?

encode possible delay

• f 1 unit by assuming

task A is length 2 the algorithm produces...

19 / 19

Lagrangian Duality → Divide & Conquer

• use duality to decompose P in smaller problems

Duality D c1 H1 A1 λ1 λ2 λI . . . xD

1

xD

2

xD

I

xD P1 P2 PI b A1 H1 H2 HI A2 AI . . . . . . . . . . . . . . . . . . c1 c2 cI . . . . . . . . . b dI . . . . . . . . . . . . c2 H2 A2 cI HI AI d1 d2 dI d1 d2

• but...

xD is generally suboptimal or even infeasible in the MIP case

• solutions may violate the relaxed constraints

i Hixi ≤ b

19 / 19

Robust PWM

ROBUST OPEN LOOP/CLOSED LOOP

19 / 19

Linear Robust MPC

Finite horizon, optimal control problem formulation min N−1

• k=0

(xk − xref )TP(xk − xref ) + uT

k Quk

• s.t.

x0 = x0, xk+1 = Axk + Buk xk ∈ Xk, uk ∈ Uk

• ideally, obtain optimal control policy using e.g. DP

◮ ”optimal decision is taken at each stage in the horizon” ◮ uk = π(xk) ◮ often intractable

• conservative approximation: open-loop control policy

◮ ”decide now for the entire horizon” (plan can’t be modifed) ◮ uk = vk ◮ poor performance, infeasibility issues

• middle ground: affine recourse

◮ ”decisions are affinely adjusted once disturb. are meas.” ◮ 19 / 19

Linear Robust MPC

Finite horizon, robust optimal control problem formulation min E N−1

• k=0

(xk − xref )TP(xk − xref ) + uT

k Quk

• s.t.

x0 = x0, xk+1 = Axk + Buk + Gwk, xk ∈ Xk, uk ∈ Uk    ∀w ∈ W

• ideally, obtain optimal control policy using e.g. DP

◮ ”optimal decision is taken at each stage in the horizon” ◮ uk = π(xk) ◮ often intractable

• conservative approximation: open-loop control policy

◮ ”decide now for the entire horizon” (plan can’t be modifed) ◮ uk = vk ◮ poor performance, infeasibility issues

• middle ground: affine recourse

◮ ”decisions are affinely adjusted once disturb. are meas.” ◮ 19 / 19

Linear Robust MPC

Finite horizon, optimal control problem formulation min E

• (x − xref)TP(x − xref) + uTQu
• s.t.

x = Ax0 + Bu + Gw, Exx + Euu ≤ e

• ∀w ∈ W
• ideally, obtain optimal control policy using e.g. DP

◮ ”optimal decision is taken at each stage in the horizon” ◮ uk = π(xk) ◮ often intractable

• conservative approximation: open-loop control policy

◮ ”decide now for the entire horizon” (plan can’t be modifed) ◮ uk = vk ◮ poor performance, infeasibility issues

• middle ground: affine recourse

◮ ”decisions are affinely adjusted once disturb. are meas.” ◮ 19 / 19

Linear Robust MPC

Finite horizon, optimal control problem formulation min E

• (x − xref)TP(x − xref) + uTQu
• s.t.

x = Ax0 + Bu + Gw, Exx + Euu ≤ e

• ∀w ∈ W
• ideally, obtain optimal control policy using e.g. DP

◮ ”optimal decision is taken at each stage in the horizon” ◮ uk = π(xk) ◮ often intractable

• conservative approximation: open-loop control policy

◮ ”decide now for the entire horizon” (plan can’t be modifed) ◮ uk = vk ◮ poor performance, infeasibility issues

• middle ground: affine recourse

◮ ”decisions are affinely adjusted once disturb. are meas.” ◮ uk = vk + k

j=0 Kkjxj

19 / 19

Linear Robust MPC

Finite horizon, optimal control problem formulation min E

• (x − xref)TP(x − xref) + uTQu
• s.t.

x = Ax0 + Bu + Gw, Exx + Euu ≤ e

• ∀w ∈ W
• ideally, obtain optimal control policy using e.g. DP

◮ ”optimal decision is taken at each stage in the horizon” ◮ uk = π(xk) ◮ often intractable

• conservative approximation: open-loop control policy

◮ ”decide now for the entire horizon” (plan can’t be modifed) ◮ uk = vk ◮ poor performance, infeasibility issues

• middle ground: affine recourse

◮ ”decisions are affinely adjusted once disturb. are meas.” ◮ uk = vk + k

j=0 Mkjwj

19 / 19

Explicit Robust Counterpart

• can show that the robust counterpart is equivalent to

minv,M,Λ E

• (x − xref)TP(x − xref) + uTQu
• s.t.

ExAx0 + (ExB + Eu)v − e ≤ −ΛT · h ΛT · S = ExG + ExBM + EuM Λ ≥ 0

(R-MPCavg) where wk is bounded by S · wk ≤ h (Goulart 2006)

• finite convex QP

19 / 19

Extension to Hybrid Systems

19 / 19

Robust Control of Hybrid Systems

• hybrid system dynamics (MLD)

xk+1 = Axk + Buk + B2δk + B3zk + Gwk Exxk + Euuk + Eδδk + Ezzk ≤ ek δk ∈ {0, 1}nδ, zk ∈ Rnz

• δk, zk characterize hybrid behavior,

◮ in the dynamics, e.g. switching between modes ◮ in the constr., e.g. logic conditions on the inputs

• wish to obtain a solution (with some recourse) to

min

x,u,δ,z

E

• (x − xref)TP(x − xref) + uTQu
• s.t.

x = Ax0 + Bu + B2δ δ δ + B3z + Gw, Exx + Euu + Eδ

δ δδ

δ δ + Ezz ≤ e , ∀w ∈ W δ δ δ ∈ {0, 1}N·nδ (RHOCP)

19 / 19

Affine recourse on continuous inputs - R-MPChyb

Main Idea

Proposed idea:

• we need recourse, but affine functions cannot easily provide binary

inputs

• split the inputs into continuous inputs u and binary inputs d

xk+1 = A · xk + Bcontuk + Bbindk + B2 · δk + B3zk + Gwk

• introduce affine recourse on the continuous inputs

u := M · w + v Assumption on G:

• disturbances only affect the continuous dynamics

19 / 19

Affine recourse on continuous inputs - R-MPChyb

• can show that the robust counterpart is

min

v,M,d,δ δ δ,z,Λ

ˆ f + trace(D · Cw) s.t. ˆ e ≤ −ΛTh , ΛTS = ExBuM + ExG + EuM Λ ≥ 0 element-wise δ δ δ ∈ {0, 1}N·nδ, d ∈ {0, 1}N·nd (R-MPChyb)

• with

ˆ x . = Ax0 + Buv + Bdd + B2δ δ δ + B3z ˆ e . = Exˆ x + Euv + Edd + Eδ

δ δδ

δ δ + Ezz − e (variables in case of zero disturbance)

• and appropriate D

19 / 19

Example: DC-DC Buck Converter

Linear Model

19 / 19

Plant: the Buck Converter (BC) (1/2)

RL RC id(t) L C Rout

Vin

iL(t) VC(t) Vo(t)

Figure: DC-DC buck converter circuit

• BC regulates input voltage Vin down to desired Vo,ref
• operated by switch (controlled input)
• disturbances: |id| ≤ 0.5 · iL,ref
• state constraints: iL ≤ 2 · iL,ref
• sampling frequency: 10 kHz

19 / 19

Averaged Model (standard method)

0.2 0.4 0.6 0.8 1 binary input control signal (A) AVERAGED cycle

• replace δ(t) ∈ {0, 1} by duty cycle d(t) ∈ [0, 1]
• average dynamics ”when off” and ”when on” weighted by d(t), obtain

˙ x = Ax + Bu + Gw u ∈ [0, 1]

• linear model

19 / 19

Performance of R-MPCavg

• 1

1 2 3

iL /iL

*

VC/VC

*

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.5 1 1.5

control signal binary input

time (ms)

19 / 19

R-MPCavg – Second Experiment

1 2 3 4 5 6 0.8 1 1.2 1.4 1.6 1.8 time (ms)

V

C,ref

constr. V

C/33V (avg)

19 / 19

Hybrid Model

Hybrid Model

19 / 19

Hybrid Model [2] (1/2)

• Closer approximation of

the switch

• Divide time into cycles with

M samples each

• New binary inputs:

switch on δ+

k = 1,

switch off δ−

k = 1

• New continuous input: uc,k
• New auxiliary state: sk

(integration of δ+ − δ−)

• switch position given by

sk + uc,k

0.2 0.4 0.6 0.8 1 binary input control signal 0.2 0.4 0.6 0.8 1 time (A) AVERAGED (B) HYBRID cycle

Figure: PWM in the hybrid model

19 / 19

Hybrid Model [2] (2/2)

• New system equations:

iL

VC s

• k+1 =

A11

A12 B1 A21 A22 B2 1

• ·

iL

VC s

• k +

B1

B2 1 −1

• ·

uc

δ+ δ−

• k +

G1

G2

• wk
• New constraints:

δ+

k , δ− k ∈ {0, 1}

binary inputs 0 ≤ sk ≤ 1 binary state 0 ≤ sk + uc,k ≤ 1 limited input −δ−

k ≤ uc,k ≤ δ+ k

switching time M−1

i=0 δ+ k+i ≤ 1

M−1

i=0 δ− k+i ≤ 1

• switching constraints
• together with iL,k ≤ 2 · iL,ref they form

xk+1 = A · xk + Bcontuk + Bbindk + B2 · δk + B3zk + Gwk Exxk + Euuk + Eddk + Eδδk + Ezzk ≤ ek

19 / 19

Performance & Simulation results – R-MPChyb

• 1

1 2 3 iL/iL

*

V

C/V C *

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.5 1 1.5 time (ms)

control signal binary input range of recourse

19 / 19

R-MPChyb – Second Experiment

1 2 3 4 5 6 0.8 1 1.2 1.4 1.6 1.8 time (ms)

V

C/33V

V

C,ref

constr. (hyb) V

C/33V (avg)

Controller RMS deviation of Vo(t) RMS deviation of Vo(t) 1st experiment 2nd experiment MPCavg 2.09V 15.2[V ] MPChyb 0.88V 9.74[V ] MPCopen loop

hyb

0.89V Infeasibility encountered

19 / 19