An Extended Mathematical Programming Framework - Semantic Scholar
An Extended Mathematical Programming Framework Michael C. Ferris University of Wisconsin, Madison
Computing with Uncertainty, IMA, October 21, 2010
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
1 / 35
Optimal Power Flow OPF(α):
miny s.t.
energy dispatch cost (y , α) conservation of power flow at nodes Kirchoff’s voltage law, and simple bound constraints
α are (given) price bids, parametric optimization
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
2 / 35
Optimal Power Flow OPF(α):
miny s.t.
energy dispatch cost (y , α) conservation of power flow at nodes Kirchoff’s voltage law, and simple bound constraints
α are (given) price bids, parametric optimization Leader(¯ α−i ):
maxαi ,y ,λ s.t.
firm i’s profit (αi , y , λ) 0 ≤ αi ≤ α ˆi y solves OPF(αi , α ¯ −i )
Note that objective involves multiplier from OPF problem
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
2 / 35
Optimal Power Flow OPF(α):
miny s.t.
energy dispatch cost (y , α) conservation of power flow at nodes Kirchoff’s voltage law, and simple bound constraints
α are (given) price bids, parametric optimization Leader(¯ α−i ):
maxαi ,y ,λ s.t.
firm i’s profit (αi , y , λ) 0 ≤ αi ≤ α ˆi y solves OPF(αi , α ¯ −i )
Note that objective involves multiplier from OPF problem Leader(¯ α−i ):
maxαi ,y ,λ s.t.
firm i’s profit (y , λ, α) 0 ≤ αi ≤ α ˆi (y , λ) solves KKT(OPF(αi , α ¯ −i ))
This is an example of an MPCC since KKT form complementarity constraints Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
2 / 35
Multi-player EPEC and security constraints (α ¯1, α ¯2, . . . , α ¯ m ) is an equilibrium if α ¯ i solves Leader(α ¯ −i ), ∀i
(Nonlinear) Nash Equilibrium where each player solves an MPCC
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
3 / 35
Multi-player EPEC and security constraints (α ¯1, α ¯2, . . . , α ¯ m ) is an equilibrium if α ¯ i solves Leader(α ¯ −i ), ∀i
(Nonlinear) Nash Equilibrium where each player solves an MPCC MPCC is hard (lacks a constraint qualification) Nash Equilibrium is PPAD-complete (Chen et al, Papadimitriou et al) In practice, also require contingency (scenario) constraints imposed in the OPF problem Model involves: complementarity, hierarchy, interacting agents, risk measures How to convey and exploit such model structure Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
3 / 35
Complementarity Problems via Graphs
T = NR+ = (R+ × {0})
S
({0} × R− )
−y ∈ T (λ) ⇐⇒ (λ, −y ) ∈ T ⇐⇒ 0 ≤ λ ⊥ y ≥ 0 By approximating (smoothing) graph can generate interior point algorithms for example y λ = , y , λ > 0 −F (x) ∈ NRn+ (x) ⇐⇒ (x, −F (x)) ∈ T n ⇐⇒ 0 ≤ x ⊥ F (x) ≥ 0 Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
4 / 35
Complementarity Systems (DVI) dx dt (t)
= f (x(t), λ(t))
y (t) = h(x(t), λ(t)) 0 ≤ y (t) ⊥ λ(t) ≥ 0
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
5 / 35
Complementarity Systems (DVI) dx dt (t)
= f (x(t), λ(t))
y (t) = h(x(t), λ(t)) 0 ≤ y (t) ⊥ λ(t) ≥ 0
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
5 / 35
Complementarity Systems (DVI) dx dt (t)
= f (x(t), λ(t))
y (t) = h(x(t), λ(t)) (λ(t), −y (t)) ∈ T
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
6 / 35
smetsyS ytiratnemelpmo
Operators and Graphs (X = [−1, 1])
xi = −1, −Fi (x) ≤ 0 or xi ∈ (−1, 1), −Fi (x) = 0 or xi = 1, −Fi (x) ≥ 0
T (λ)
Ferris (Univ. Wisconsin)
T −1 (y ) EMP
PT (y ) is the projection of y onto [`, u] Ferris (Univ. Wisconsin)
EMP
(I + T )−1 (y ) = PT (y ) Ferris Nonsmooth (Univ. Wisconsin) School, June 2010
IMA, October 2010
7 / 35
2
Generalized Equations Suppose T is a maximal monotone operator 0 ∈ F (z) + T (z)
(GE )
Define PT = (I + T )−1 If T is polyhedral (graph of T is a finite union of convex polyhedral sets) then PT is piecewise affine (continous, single-valued, non-expansive) 0 ∈ F (z) + T (z)
⇐⇒
z ∈ F (z) + I(z) + T (z)
⇐⇒
z − F (z) ∈ (I + T )(z) ⇐⇒ PT (z − F (z)) = z
Use in fixed point iterations (cf projected gradient methods)
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
8 / 35
Splitting Methods Suppose T is a maximal monotone operator 0 ∈ F (z) + T (z)
(GE )
Can devise Newton methods (e.g. SQP) that treat F via calculus and T via convex analysis Alternatively, can split F (z) = A(z) + B(z) (and possibly T also) so we solve solve (GE) by solving a sequence of problems involving just T1 (z) = A(z) and T2 (z) = B(z) + T (z) where each of these is “simpler” Forward-Backward splitting: z k+1 = (I + ck T2 )−1 (I − ck T1 ) z k , Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
9 / 35
Normal Map Suppose T is a maximal monotone operator 0 ∈ F (z) + T (z)
(GE )
Define PT = (I + T )−1 0 ∈ F (z) + T (z)
⇐⇒
z ∈ F (z) + I(z) + T (z)
⇐⇒
z − F (z) = y and y ∈ (I + T )(z)
⇐⇒
z − F (z) = y and PT (y ) = z
⇐⇒
PT (y ) − F (PT (y )) = y
⇐⇒
0 = F (PT (y )) + y − PT (y )
This is the so-called Normal Map Equation
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
10 / 35
Normal manifold = {Fi + NFi }
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
11 / 35
C = {z|Bz ≥ b}, F (z) = Mz + q
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
12 / 35
Cao/Ferris Path (Eaves)
Cao/Ferris Path
Start in cell that has interior (face is an extreme point) Move towards a zero of affine map in cell Update direction when hit boundary (pivot) Solves or determines infeasible if M is copositive-plus on rec(C )
Solves 2-person bimatrix games, 3-person games too, but these are nonlinear But algorithm has exponential complexity (von Stengel et al)
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
• Sta inte ext • Mov affi • Upd hit • Solv infe cop • Nai
13 / 35
Extensions and Computational Results Embed AVI solver in a Newton Method - each Newton step solves an AVI Compare performance of PathAVI with PATH on equivalent LCP PATH the most widely used code for solving MCP AVIs constructed to have solution with Mn×n symmetric indefinite
Size (m,n) (180, 60) (360, 120)
PathAVI Resid Iter 3 × 10−14 193 3 × 10−14 1516
PATH Resid 0.9 2.0
Iter 10176 10594
2 - 10x speedup in Matlab using sparse LU instead of QR 2 - 10x speedup in C using sparse LU updates
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
14 / 35
But who cares?
Why aren’t you using my *********** algorithm? (Michael Ferris, Boulder, CO, 1994)
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
15 / 35
But who cares?
Why aren’t you using my *********** algorithm? (Michael Ferris, Boulder, CO, 1994) Show me on a problem like mine Must run on defaults Must deal graciously with poorly specified cases Must be usable from my environment (Matlab, R, GAMS, ...) Must be able to model my problem easily EMP provides annotations to an existing optimization model that convey new model structures to a solver
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
15 / 35
EMP(i): MPCC: complementarity constraints min x,s
s.t.
f (x, s) g (x, s) ≤ 0, 0 ≤ s ⊥ h(x, s) ≥ 0
g , h model “engineering” expertise: finite elements, etc ⊥ models complementarity, disjunctions Complementarity “⊥” constraints available in AIMMS, AMPL and GAMS
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
16 / 35
EMP(i): MPCC: complementarity constraints min x,s
s.t.
f (x, s) g (x, s) ≤ 0, 0 ≤ s ⊥ h(x, s) ≥ 0
g , h model “engineering” expertise: finite elements, etc ⊥ models complementarity, disjunctions Complementarity “⊥” constraints available in AIMMS, AMPL and GAMS NLPEC: use the convert tool to automatically reformulate as a parameteric sequence of NLP’s Solution by repeated use of standard NLP software I
Problems solvable, local solutions, hard
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
16 / 35
EMP(ii): Hierarchical models Bilevel programs: min
f (x ∗ , y ∗ )
s.t.
g (x ∗ , y ∗ ) ≤ 0, y ∗ solves min v (x ∗ , y ) s.t. h(x ∗ , y ) ≤ 0
x ∗ ,y ∗
y
model bilev /deff,defg,defv,defh/; empinfo: bilevel min v y defv defh EMP tool automatically creates the MPCC min
x ∗ ,y ∗ ,λ
s.t.
Ferris (Univ. Wisconsin)
f (x ∗ , y ∗ ) g (x ∗ , y ∗ ) ≤ 0, 0 ≤ ∇v (x ∗ , y ∗ ) + λT ∇h(x ∗ , y ∗ ) ⊥ y ∗ ≥ 0 0 ≤ −h(x ∗ , y ∗ ) ⊥λ≥0 EMP
IMA, October 2010
17 / 35
Large scale example: bioreactor (with Niebel)
Challenge Formulating an optimization problem that allows the estimation of the dynamic changes in intracellular fluxes based on measured external bioreactor concentrations.
Approach Using existing constraint-based stoichiometric models of the cellular metabolism to formulate a bilevel dynamic optimization problem.
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
18 / 35
Dynamic optimization Approach: The different timescales of the metabolism (fast) and the reactor growth (slow), allows to assume steady-state for the metabolism. minimize / maximize Objective (eg. parameter fitting) s. t.
bioreactor dynamics constraints on exchange fluxes maximize growth rate s. t. stoichiometric constraints flux constraints
Different mathematical programming techniques are used to transform the problem to a nonlinear program. The differential equations are transformed into nonlinear constraints using collocation methods. Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
19 / 35
The overall scheme!
Collection of algebraic equations Form a bilevel program via emp EMP tool automatically creates the MPCC (model transformation) NLPEC tool automatically creates (a series of) NLP models (model transformation) GAMS automatically rewrites NLP models for global solution via BARON (model transformation) Is this global? What’s the hitch?
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
20 / 35
The overall scheme!
Collection of algebraic equations Form a bilevel program via emp EMP tool automatically creates the MPCC (model transformation) NLPEC tool automatically creates (a series of) NLP models (model transformation) GAMS automatically rewrites NLP models for global solution via BARON (model transformation) Is this global? What’s the hitch? Note that heirarchical structure is available to solvers for analysis or utilization
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
20 / 35
Variational inequalities Find z ∈ X such that 0 ∈ F (z) + NX (z)
Many applications where F is not the derivative of some f model vi / F, g /; empinfo: vifunc F z Convert problem into complementarity problem by introducing multipliers on representation of X Can now do MPEC (as opposed to MPCC)! Projection algorithms, robustness (evaluate F only at points in X )
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
21 / 35
EMP(iii): Embedded models Model has the format: Agent o: min f (x, y ) x
s.t.
g (x, y ) ≤ 0 (⊥ λ ≥ 0)
Agent v: H(x, y , λ) = 0
(⊥ y free)
Difficult to implement correctly (multiple optimization models) Can do automatically - simply annotate equations empinfo: equilibrium min f x defg vifunc H y dualvar λ defg EMP tool automatically creates an MCP ∇x f (x, y ) + λT ∇g (x, y ) = 0 0 ≤ −g (x, y ) ⊥ λ ≥ 0 H(x, y , λ) = 0 Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
22 / 35
World Bank Project (Uruguay Round)
Application: Uruguay Round
24 regions, 22 commodities
Nonlinear complementarity problem I Size:Bank 2200Project x 2200 with • World Harrison and Rutherford Short term gains $53 billion p.a. • 24 regions, 22 commodities I Much smaller than previous I
– 2200 x 2200 (nonlinear)
literature • Short term gains $53 billion p.a.
Long term gains $188 billion p.a. – Much smaller than previous
Number literatureof less developed countries loose$188 in short termp.a. • Long term gains billion I
– Numberconclusions of less developed Unpopular - forced countries loose in short term concessions by World Bank • Unpopular conclusions – forced
concessions by World Bank not Region/commodity structure apparent to solver
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
23 / 35
Nash Equilibria Nash Games: x ∗ is a Nash Equilibrium if ∗ xi∗ ∈ arg min `i (xi , x−i , q), ∀i ∈ I xi ∈Xi
x−i are the decisions of other players. Quantities q given exogenously, or via complementarity: 0 ≤ H(x, q) ⊥ q ≥ 0
empinfo: equilibrium min loss(i) x(i) cons(i) vifunc H q Applications: Discrete-Time Finite-State Stochastic Games. Specifically, the Ericson & Pakes (1995) model of dynamic competition in an oligopolistic industry. Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
24 / 35
Key point: models generated correctly solve quickly Here S is mesh spacing parameter S Var rows non-zero dense(%) Steps RT (m:s) 20 2400 2568 31536 0.48 5 0 : 03 50 15000 15408 195816 0.08 5 0 : 19 100 60000 60808 781616 0.02 5 1 : 16 200 240000 241608 3123216 0.01 5 5 : 12 Convergence for S = 200 (with new basis extensions in PATH) Iteration 0 1 2 3 4 5 Ferris (Univ. Wisconsin)
Residual 1.56(+4) 1.06(+1) 1.34 2.04(−2) 1.74(−5) 2.97(−11) EMP
IMA, October 2010
25 / 35
General Equilibrium models (C ) : max Uk (xk ) s.t. p T xk ≤ ik (y , p) xk ∈Xk X (I ) :ik (y , p) = p T ωk + αkj p T gj (yj ) j T
(P) : max p gj (yj ) yj ∈Yj X X X X (M) : max p T xk − ωk − gj (yj ) s.t. pl = 1 p≥0
Ferris (Univ. Wisconsin)
k
j
k
EMP
l
IMA, October 2010
26 / 35
General Equilibrium models (C ) : max Uk (xk ) s.t. p T xk ≤ ik (y , p) xk ∈Xk X (I ) :ik (y , p) = p T ωk + αkj p T gj (yj ) j T
(P) : max p gj (yj ) yj ∈Yj X X X X (M) : max p T xk − ωk − gj (yj ) s.t. pl = 1 p≥0
k
j
k
l
Can reformulate as embedded problem (Ermoliev et al): X tk log Uk (xk ) max x∈X ,y ∈Y βk k X X X s.t. xk ≤ ωk + gj (yj ) k
j
k
tk = ik (y , p) where p is multiplier on NLP constraint Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
26 / 35
Competing agent models (with Wets)
Competing agents (consumers, or generators in energy market) Each agent maximizes objective independently (utility) Market prices are function of all agents activities Additional twist: model must “hedge” against uncertainty Facilitated by allowing contracts bought now, for goods delivered later Conceptually allows to transfer goods from one period to another (provides wealth retention or pricing of ancilliary services in energy market)
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
27 / 35
The model details: c.f. Brown, Demarzo, Eaves Each agent maximizes: ! uh = −
X
κ−
πs
s
Y
αh,l ch,s,l
l
Time 0: X
dh,0,l = ch,0,l − eh,0,l ,
p0,l dh,0,l +
X
l
qk zh,k ≤ 0
k
Time 1: dh,s,l = ch,s,l − eh,s,l −
X
Ds,l,k ∗ zh,k ,
X
k
ps,l dh,s,l ≤ 0
l
Additional constraints (complementarity) outside of control of agents: 0≤−
X
zh,k ⊥ qk ≥ 0
h
0≤− Ferris (Univ. Wisconsin)
X h
dh,s,l ⊥ ps,l ≥ 0 EMP
IMA, October 2010
28 / 35
EMP(iv): Stochastic programming and risk measures SP: min c > x + E[d > y ] s.t. Ax = b T (ω)x + W (ω)y (ω) ≥ h(ω), x ≥ 0,
y (ω) ≥ 0,
for all ω ∈ Ω,
for all ω ∈ Ω.
Annotations are slightly more involved but straightforward: Need to describe probability distribution Define (multi-stage) structure (what variables and constraints belong to each stage) Define random parameters and process to generate scenarios Can also define risk measures on variables
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
29 / 35
EMP(iv): Stochastic programming and risk measures SP: min c > x + R[d > y ] s.t. Ax = b T (ω)x + W (ω)y (ω) ≥ h(ω), x ≥ 0,
y (ω) ≥ 0,
for all ω ∈ Ω,
for all ω ∈ Ω.
Annotations are slightly more involved but straightforward: Need to describe probability distribution Define (multi-stage) structure (what variables and constraints belong to each stage) Define random parameters and process to generate scenarios Can also define risk measures on variables Automatic reformulation (deterministic equivalent), solvers such as DECIS, etc. Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
30 / 35
EMP(v): Extended nonlinear programs min f0 (x)+θ(f1 (x), . . . , fm (x))
x∈X
Examples of different θ
least squares, absolute value, Huber function Solution reformulations are very different Huber function used in robust statistics.
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
31 / 35
More general θ functions
In general any piecewise linear penalty function can be used: (different upside/downside costs). General form: θ(u) = sup {y T u − k(y )} y ∈Y
θ nonsmooth due to the max term; θ separable in example. θ is always convex. Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
32 / 35
First order conditions
Solution via reformulation. One way: 0 ∈ ∇x L(x, y ) + NX (x) 0 ∈ −∇y L(x, y ) + NY (y ) NX (x) is the normal cone to the closed convex set X at x. Automatically creates an MCP (or a VI) Already available! To do: extend X and Y beyond simple bound sets.
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
33 / 35
Alternative Reformulations Convert does symbolic/numeric reformulations. Alternative NLP formulations also possible. 1 k(y ) = y 0 Qy , X = {x : Rx ≤ r } , Y = y : S 0 y ≤ s 2 Defining Q = DJ −1 D 0 , F (x) = (f1 (x), . . . , fm (x)) min f0 (x) + s 0 z + 12 wJw s.t. Rx ≤ r , z ≥ 0, F (x) − Sz − Dw = 0 Can set up better (solver) specific formulation.
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
34 / 35
Conclusions Modern optimization within applications requires multiple model formats, computational tools and sophisticated solvers EMP model type is clear and extensible, additional structure available to solver Extended Mathematical Programming available within the GAMS modeling system Able to pass additional (structure) information to solvers Embedded optimization models automatically reformulated for appropriate solution engine Exploit structure in solvers Extend application usage further
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
35 / 35
Computing with Uncertainty, IMA, October 21, 2010
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
1 / 35
Optimal Power Flow OPF(α):
miny s.t.
energy dispatch cost (y , α) conservation of power flow at nodes Kirchoff’s voltage law, and simple bound constraints
α are (given) price bids, parametric optimization
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
2 / 35
Optimal Power Flow OPF(α):
miny s.t.
energy dispatch cost (y , α) conservation of power flow at nodes Kirchoff’s voltage law, and simple bound constraints
α are (given) price bids, parametric optimization Leader(¯ α−i ):
maxαi ,y ,λ s.t.
firm i’s profit (αi , y , λ) 0 ≤ αi ≤ α ˆi y solves OPF(αi , α ¯ −i )
Note that objective involves multiplier from OPF problem
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
2 / 35
Optimal Power Flow OPF(α):
miny s.t.
energy dispatch cost (y , α) conservation of power flow at nodes Kirchoff’s voltage law, and simple bound constraints
α are (given) price bids, parametric optimization Leader(¯ α−i ):
maxαi ,y ,λ s.t.
firm i’s profit (αi , y , λ) 0 ≤ αi ≤ α ˆi y solves OPF(αi , α ¯ −i )
Note that objective involves multiplier from OPF problem Leader(¯ α−i ):
maxαi ,y ,λ s.t.
firm i’s profit (y , λ, α) 0 ≤ αi ≤ α ˆi (y , λ) solves KKT(OPF(αi , α ¯ −i ))
This is an example of an MPCC since KKT form complementarity constraints Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
2 / 35
Multi-player EPEC and security constraints (α ¯1, α ¯2, . . . , α ¯ m ) is an equilibrium if α ¯ i solves Leader(α ¯ −i ), ∀i
(Nonlinear) Nash Equilibrium where each player solves an MPCC
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
3 / 35
Multi-player EPEC and security constraints (α ¯1, α ¯2, . . . , α ¯ m ) is an equilibrium if α ¯ i solves Leader(α ¯ −i ), ∀i
(Nonlinear) Nash Equilibrium where each player solves an MPCC MPCC is hard (lacks a constraint qualification) Nash Equilibrium is PPAD-complete (Chen et al, Papadimitriou et al) In practice, also require contingency (scenario) constraints imposed in the OPF problem Model involves: complementarity, hierarchy, interacting agents, risk measures How to convey and exploit such model structure Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
3 / 35
Complementarity Problems via Graphs
T = NR+ = (R+ × {0})
S
({0} × R− )
−y ∈ T (λ) ⇐⇒ (λ, −y ) ∈ T ⇐⇒ 0 ≤ λ ⊥ y ≥ 0 By approximating (smoothing) graph can generate interior point algorithms for example y λ = , y , λ > 0 −F (x) ∈ NRn+ (x) ⇐⇒ (x, −F (x)) ∈ T n ⇐⇒ 0 ≤ x ⊥ F (x) ≥ 0 Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
4 / 35
Complementarity Systems (DVI) dx dt (t)
= f (x(t), λ(t))
y (t) = h(x(t), λ(t)) 0 ≤ y (t) ⊥ λ(t) ≥ 0
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
5 / 35
Complementarity Systems (DVI) dx dt (t)
= f (x(t), λ(t))
y (t) = h(x(t), λ(t)) 0 ≤ y (t) ⊥ λ(t) ≥ 0
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
5 / 35
Complementarity Systems (DVI) dx dt (t)
= f (x(t), λ(t))
y (t) = h(x(t), λ(t)) (λ(t), −y (t)) ∈ T
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
6 / 35
smetsyS ytiratnemelpmo
Operators and Graphs (X = [−1, 1])
xi = −1, −Fi (x) ≤ 0 or xi ∈ (−1, 1), −Fi (x) = 0 or xi = 1, −Fi (x) ≥ 0
T (λ)
Ferris (Univ. Wisconsin)
T −1 (y ) EMP
PT (y ) is the projection of y onto [`, u] Ferris (Univ. Wisconsin)
EMP
(I + T )−1 (y ) = PT (y ) Ferris Nonsmooth (Univ. Wisconsin) School, June 2010
IMA, October 2010
7 / 35
2
Generalized Equations Suppose T is a maximal monotone operator 0 ∈ F (z) + T (z)
(GE )
Define PT = (I + T )−1 If T is polyhedral (graph of T is a finite union of convex polyhedral sets) then PT is piecewise affine (continous, single-valued, non-expansive) 0 ∈ F (z) + T (z)
⇐⇒
z ∈ F (z) + I(z) + T (z)
⇐⇒
z − F (z) ∈ (I + T )(z) ⇐⇒ PT (z − F (z)) = z
Use in fixed point iterations (cf projected gradient methods)
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
8 / 35
Splitting Methods Suppose T is a maximal monotone operator 0 ∈ F (z) + T (z)
(GE )
Can devise Newton methods (e.g. SQP) that treat F via calculus and T via convex analysis Alternatively, can split F (z) = A(z) + B(z) (and possibly T also) so we solve solve (GE) by solving a sequence of problems involving just T1 (z) = A(z) and T2 (z) = B(z) + T (z) where each of these is “simpler” Forward-Backward splitting: z k+1 = (I + ck T2 )−1 (I − ck T1 ) z k , Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
9 / 35
Normal Map Suppose T is a maximal monotone operator 0 ∈ F (z) + T (z)
(GE )
Define PT = (I + T )−1 0 ∈ F (z) + T (z)
⇐⇒
z ∈ F (z) + I(z) + T (z)
⇐⇒
z − F (z) = y and y ∈ (I + T )(z)
⇐⇒
z − F (z) = y and PT (y ) = z
⇐⇒
PT (y ) − F (PT (y )) = y
⇐⇒
0 = F (PT (y )) + y − PT (y )
This is the so-called Normal Map Equation
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
10 / 35
Normal manifold = {Fi + NFi }
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
11 / 35
C = {z|Bz ≥ b}, F (z) = Mz + q
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
12 / 35
Cao/Ferris Path (Eaves)
Cao/Ferris Path
Start in cell that has interior (face is an extreme point) Move towards a zero of affine map in cell Update direction when hit boundary (pivot) Solves or determines infeasible if M is copositive-plus on rec(C )
Solves 2-person bimatrix games, 3-person games too, but these are nonlinear But algorithm has exponential complexity (von Stengel et al)
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
• Sta inte ext • Mov affi • Upd hit • Solv infe cop • Nai
13 / 35
Extensions and Computational Results Embed AVI solver in a Newton Method - each Newton step solves an AVI Compare performance of PathAVI with PATH on equivalent LCP PATH the most widely used code for solving MCP AVIs constructed to have solution with Mn×n symmetric indefinite
Size (m,n) (180, 60) (360, 120)
PathAVI Resid Iter 3 × 10−14 193 3 × 10−14 1516
PATH Resid 0.9 2.0
Iter 10176 10594
2 - 10x speedup in Matlab using sparse LU instead of QR 2 - 10x speedup in C using sparse LU updates
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
14 / 35
But who cares?
Why aren’t you using my *********** algorithm? (Michael Ferris, Boulder, CO, 1994)
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
15 / 35
But who cares?
Why aren’t you using my *********** algorithm? (Michael Ferris, Boulder, CO, 1994) Show me on a problem like mine Must run on defaults Must deal graciously with poorly specified cases Must be usable from my environment (Matlab, R, GAMS, ...) Must be able to model my problem easily EMP provides annotations to an existing optimization model that convey new model structures to a solver
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
15 / 35
EMP(i): MPCC: complementarity constraints min x,s
s.t.
f (x, s) g (x, s) ≤ 0, 0 ≤ s ⊥ h(x, s) ≥ 0
g , h model “engineering” expertise: finite elements, etc ⊥ models complementarity, disjunctions Complementarity “⊥” constraints available in AIMMS, AMPL and GAMS
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
16 / 35
EMP(i): MPCC: complementarity constraints min x,s
s.t.
f (x, s) g (x, s) ≤ 0, 0 ≤ s ⊥ h(x, s) ≥ 0
g , h model “engineering” expertise: finite elements, etc ⊥ models complementarity, disjunctions Complementarity “⊥” constraints available in AIMMS, AMPL and GAMS NLPEC: use the convert tool to automatically reformulate as a parameteric sequence of NLP’s Solution by repeated use of standard NLP software I
Problems solvable, local solutions, hard
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
16 / 35
EMP(ii): Hierarchical models Bilevel programs: min
f (x ∗ , y ∗ )
s.t.
g (x ∗ , y ∗ ) ≤ 0, y ∗ solves min v (x ∗ , y ) s.t. h(x ∗ , y ) ≤ 0
x ∗ ,y ∗
y
model bilev /deff,defg,defv,defh/; empinfo: bilevel min v y defv defh EMP tool automatically creates the MPCC min
x ∗ ,y ∗ ,λ
s.t.
Ferris (Univ. Wisconsin)
f (x ∗ , y ∗ ) g (x ∗ , y ∗ ) ≤ 0, 0 ≤ ∇v (x ∗ , y ∗ ) + λT ∇h(x ∗ , y ∗ ) ⊥ y ∗ ≥ 0 0 ≤ −h(x ∗ , y ∗ ) ⊥λ≥0 EMP
IMA, October 2010
17 / 35
Large scale example: bioreactor (with Niebel)
Challenge Formulating an optimization problem that allows the estimation of the dynamic changes in intracellular fluxes based on measured external bioreactor concentrations.
Approach Using existing constraint-based stoichiometric models of the cellular metabolism to formulate a bilevel dynamic optimization problem.
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
18 / 35
Dynamic optimization Approach: The different timescales of the metabolism (fast) and the reactor growth (slow), allows to assume steady-state for the metabolism. minimize / maximize Objective (eg. parameter fitting) s. t.
bioreactor dynamics constraints on exchange fluxes maximize growth rate s. t. stoichiometric constraints flux constraints
Different mathematical programming techniques are used to transform the problem to a nonlinear program. The differential equations are transformed into nonlinear constraints using collocation methods. Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
19 / 35
The overall scheme!
Collection of algebraic equations Form a bilevel program via emp EMP tool automatically creates the MPCC (model transformation) NLPEC tool automatically creates (a series of) NLP models (model transformation) GAMS automatically rewrites NLP models for global solution via BARON (model transformation) Is this global? What’s the hitch?
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
20 / 35
The overall scheme!
Collection of algebraic equations Form a bilevel program via emp EMP tool automatically creates the MPCC (model transformation) NLPEC tool automatically creates (a series of) NLP models (model transformation) GAMS automatically rewrites NLP models for global solution via BARON (model transformation) Is this global? What’s the hitch? Note that heirarchical structure is available to solvers for analysis or utilization
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
20 / 35
Variational inequalities Find z ∈ X such that 0 ∈ F (z) + NX (z)
Many applications where F is not the derivative of some f model vi / F, g /; empinfo: vifunc F z Convert problem into complementarity problem by introducing multipliers on representation of X Can now do MPEC (as opposed to MPCC)! Projection algorithms, robustness (evaluate F only at points in X )
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
21 / 35
EMP(iii): Embedded models Model has the format: Agent o: min f (x, y ) x
s.t.
g (x, y ) ≤ 0 (⊥ λ ≥ 0)
Agent v: H(x, y , λ) = 0
(⊥ y free)
Difficult to implement correctly (multiple optimization models) Can do automatically - simply annotate equations empinfo: equilibrium min f x defg vifunc H y dualvar λ defg EMP tool automatically creates an MCP ∇x f (x, y ) + λT ∇g (x, y ) = 0 0 ≤ −g (x, y ) ⊥ λ ≥ 0 H(x, y , λ) = 0 Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
22 / 35
World Bank Project (Uruguay Round)
Application: Uruguay Round
24 regions, 22 commodities
Nonlinear complementarity problem I Size:Bank 2200Project x 2200 with • World Harrison and Rutherford Short term gains $53 billion p.a. • 24 regions, 22 commodities I Much smaller than previous I
– 2200 x 2200 (nonlinear)
literature • Short term gains $53 billion p.a.
Long term gains $188 billion p.a. – Much smaller than previous
Number literatureof less developed countries loose$188 in short termp.a. • Long term gains billion I
– Numberconclusions of less developed Unpopular - forced countries loose in short term concessions by World Bank • Unpopular conclusions – forced
concessions by World Bank not Region/commodity structure apparent to solver
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
23 / 35
Nash Equilibria Nash Games: x ∗ is a Nash Equilibrium if ∗ xi∗ ∈ arg min `i (xi , x−i , q), ∀i ∈ I xi ∈Xi
x−i are the decisions of other players. Quantities q given exogenously, or via complementarity: 0 ≤ H(x, q) ⊥ q ≥ 0
empinfo: equilibrium min loss(i) x(i) cons(i) vifunc H q Applications: Discrete-Time Finite-State Stochastic Games. Specifically, the Ericson & Pakes (1995) model of dynamic competition in an oligopolistic industry. Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
24 / 35
Key point: models generated correctly solve quickly Here S is mesh spacing parameter S Var rows non-zero dense(%) Steps RT (m:s) 20 2400 2568 31536 0.48 5 0 : 03 50 15000 15408 195816 0.08 5 0 : 19 100 60000 60808 781616 0.02 5 1 : 16 200 240000 241608 3123216 0.01 5 5 : 12 Convergence for S = 200 (with new basis extensions in PATH) Iteration 0 1 2 3 4 5 Ferris (Univ. Wisconsin)
Residual 1.56(+4) 1.06(+1) 1.34 2.04(−2) 1.74(−5) 2.97(−11) EMP
IMA, October 2010
25 / 35
General Equilibrium models (C ) : max Uk (xk ) s.t. p T xk ≤ ik (y , p) xk ∈Xk X (I ) :ik (y , p) = p T ωk + αkj p T gj (yj ) j T
(P) : max p gj (yj ) yj ∈Yj X X X X (M) : max p T xk − ωk − gj (yj ) s.t. pl = 1 p≥0
Ferris (Univ. Wisconsin)
k
j
k
EMP
l
IMA, October 2010
26 / 35
General Equilibrium models (C ) : max Uk (xk ) s.t. p T xk ≤ ik (y , p) xk ∈Xk X (I ) :ik (y , p) = p T ωk + αkj p T gj (yj ) j T
(P) : max p gj (yj ) yj ∈Yj X X X X (M) : max p T xk − ωk − gj (yj ) s.t. pl = 1 p≥0
k
j
k
l
Can reformulate as embedded problem (Ermoliev et al): X tk log Uk (xk ) max x∈X ,y ∈Y βk k X X X s.t. xk ≤ ωk + gj (yj ) k
j
k
tk = ik (y , p) where p is multiplier on NLP constraint Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
26 / 35
Competing agent models (with Wets)
Competing agents (consumers, or generators in energy market) Each agent maximizes objective independently (utility) Market prices are function of all agents activities Additional twist: model must “hedge” against uncertainty Facilitated by allowing contracts bought now, for goods delivered later Conceptually allows to transfer goods from one period to another (provides wealth retention or pricing of ancilliary services in energy market)
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
27 / 35
The model details: c.f. Brown, Demarzo, Eaves Each agent maximizes: ! uh = −
X
κ−
πs
s
Y
αh,l ch,s,l
l
Time 0: X
dh,0,l = ch,0,l − eh,0,l ,
p0,l dh,0,l +
X
l
qk zh,k ≤ 0
k
Time 1: dh,s,l = ch,s,l − eh,s,l −
X
Ds,l,k ∗ zh,k ,
X
k
ps,l dh,s,l ≤ 0
l
Additional constraints (complementarity) outside of control of agents: 0≤−
X
zh,k ⊥ qk ≥ 0
h
0≤− Ferris (Univ. Wisconsin)
X h
dh,s,l ⊥ ps,l ≥ 0 EMP
IMA, October 2010
28 / 35
EMP(iv): Stochastic programming and risk measures SP: min c > x + E[d > y ] s.t. Ax = b T (ω)x + W (ω)y (ω) ≥ h(ω), x ≥ 0,
y (ω) ≥ 0,
for all ω ∈ Ω,
for all ω ∈ Ω.
Annotations are slightly more involved but straightforward: Need to describe probability distribution Define (multi-stage) structure (what variables and constraints belong to each stage) Define random parameters and process to generate scenarios Can also define risk measures on variables
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
29 / 35
EMP(iv): Stochastic programming and risk measures SP: min c > x + R[d > y ] s.t. Ax = b T (ω)x + W (ω)y (ω) ≥ h(ω), x ≥ 0,
y (ω) ≥ 0,
for all ω ∈ Ω,
for all ω ∈ Ω.
Annotations are slightly more involved but straightforward: Need to describe probability distribution Define (multi-stage) structure (what variables and constraints belong to each stage) Define random parameters and process to generate scenarios Can also define risk measures on variables Automatic reformulation (deterministic equivalent), solvers such as DECIS, etc. Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
30 / 35
EMP(v): Extended nonlinear programs min f0 (x)+θ(f1 (x), . . . , fm (x))
x∈X
Examples of different θ
least squares, absolute value, Huber function Solution reformulations are very different Huber function used in robust statistics.
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
31 / 35
More general θ functions
In general any piecewise linear penalty function can be used: (different upside/downside costs). General form: θ(u) = sup {y T u − k(y )} y ∈Y
θ nonsmooth due to the max term; θ separable in example. θ is always convex. Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
32 / 35
First order conditions
Solution via reformulation. One way: 0 ∈ ∇x L(x, y ) + NX (x) 0 ∈ −∇y L(x, y ) + NY (y ) NX (x) is the normal cone to the closed convex set X at x. Automatically creates an MCP (or a VI) Already available! To do: extend X and Y beyond simple bound sets.
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
33 / 35
Alternative Reformulations Convert does symbolic/numeric reformulations. Alternative NLP formulations also possible. 1 k(y ) = y 0 Qy , X = {x : Rx ≤ r } , Y = y : S 0 y ≤ s 2 Defining Q = DJ −1 D 0 , F (x) = (f1 (x), . . . , fm (x)) min f0 (x) + s 0 z + 12 wJw s.t. Rx ≤ r , z ≥ 0, F (x) − Sz − Dw = 0 Can set up better (solver) specific formulation.
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
34 / 35
Conclusions Modern optimization within applications requires multiple model formats, computational tools and sophisticated solvers EMP model type is clear and extensible, additional structure available to solver Extended Mathematical Programming available within the GAMS modeling system Able to pass additional (structure) information to solvers Embedded optimization models automatically reformulated for appropriate solution engine Exploit structure in solvers Extend application usage further
Ferris (Univ. Wisconsin)
EMP
IMA, October 2010
35 / 35
