Programmazione Lineare
INTERA III
Generazione di colonne
riferimento per lo studio: nel folder “materiale didattico integrativo”
nel minisito del corso, scaricare il file “Generazione_di_Colonne.pdf”
•
Quando la formulazione migliore (il cui rilassamento continuo ha un integrality gap molto ridotto)
presenta un numero ingestibile (esponenziale) di variabili.
E’ possibile gestirle dinamicamente?
•
Osserviamo che nella soluzione di base ottima di un problema di minimo le variabili fuori base
(=0) hanno costo ridotto 0,  posso non esplicitare le colonne che corrispondono a tali variabili

posso iterativamente individuare le colonne non esplicitate con costo ridotto <0 e aggiungerle
oppure - se non ne esistono - certificare l’ottimo del rilassamento continuo
•
•
Osservazione:
il procedimento di “generazione di colonne” e’ molto simile a quanto fatto per la soluzione del
rilassamento continuo di quei problemi la cui formulazione e’ caratterizzata da un numero molto
elevato di vincoli, sia che siano necessari alla correttezza del modello sia che siano utili a
rafforzare il rilassamento continuo, per rendere il politopo assocoato il piu’ possibile aderente al
politopo della formulazione ideale (quella i cui punti estremi hanno coordinate intere, e il cui
politopo associato corrisponde all’inviluppo convesso delle soluzini itere del problema).
Idea: in base alla teoria della dualita’ nella programmazione lineare, generare colonne nel primale
equivale a generare vincoli nel problema duale.
Vedremo come generare una colonna che ha costo ridotto negativo (in un primale di minimo)
equivale matematicamente a generare un vincolo violato dalla soluzione continua corrente nel
problema duale, in base al principio che primale ottimalita’  duale ammissibilita’
Funziona?
• Solo se posso trasformare la ricerca delle
colonne mancanti in un nuovo problema di
ottimizzazione che so risolvere in modo
efficiente (senza la ricerca esaustiva).
• Funzione obiettivo: minimizza costo ridotto della
nuova colonna
• Vincolo di ammissibilita’: il vettore soluzione
deve essere una colonna del problema
Bin Packing
Dati
• dei contenitori (bin) J di capacita’ Q
• Un set di oggetti I ={1,..,n} di peso {pi}
Si vuole minimizzare il numero di bin
necessari a contenere TUTTI gli oggetti
Modello matematico di PLI
Un possibile modello matematico utilizza:
• delle variabili y per l’uso dei bin
• delle variabili x per assegnare gli oggetti ai bin.
Min jJ yj
jJ xij = 1
 iI pijxij  Qyj
xij{0,1}
yj{0,1}
 iI
 jJ
 iI jJ
 jJ
Lower Bound del rilassamento continuo?
 ipi/Q 
Si tratta di un modello con integrality gap molto alto
E’ come se nel rilassamento continuo potessi giustapporre i bin l’uno all’altro senza soluzione
di continuita’ e ripartire un singolo oggetto fra + bin contigui, e considerare come il costo di
ciascun bin solo la frazione occupata
Posso migliorare il modello precedente aggiungendo esplicitamente i tagli xijyj
ma resta un modello inefficiente a causa di un elevato grado di simmetria nella soluzione:
Infatti:
Per ogni soluzione con k bin, esistono k! soluzioni equivalenti dal punto di vista della
ripartizione degli oggetti nei singoli bin, e che differiscono solo per la numerazione data ai
singoli bin
Min jJ yj
jJ xij = 1
 iI pijxij  Q
xij  yj
xij{0,1}
yj{0,1}
 iI
 jJ
 iI, jJ
 iI jJ
 jJ
esercizio
• Utilizzando Xpress, confronta il valore ottimo del
rilassamento continuo delle possibili
formulazioni alternative di PLI del problema di
Bin Packing per le istanze:
P={7,2,3,4,6,1}
Q=8
P={8,3,4,8,4} Q=14
Alternative?
Un nuovo modello….
Un nuovo punto di vista….
Enumeriamo i possibili modi di impaccare
degli oggetti in un bin P={1,2,3,4,6,7} Q=8
Q=8
3
4
7
2
1
3
2
6
Generiamo TUTTI i “pattern di impaccamento” possibili per P={1,2,3,4,6,7},
Q=8. Qui sotto sono quasi tutti rappresentati, tranne alcuni di quelli formati da
un singolo pezzo, i cosidetti “pattern elementari”, dei pezzi di lunghezza 1,2,3,4
1
4
4
3
7
1
1
1
4
6
2
2
2
1
4
4
4
3
7
2
1
3
6
2
6
Esercizio:
se ordino per dimensione crescente gli elementi nell’insieme degli oggetti,
posso usare il vettore caratteristico per identificare ciascun pattern.
Infatti l’ordine in cui i singoli pezzi comapiono nel pattern NON e’ rilevante,
ma conta solo se ciascun elemento Iesimo e’ o meno presente (valore 0 o 1
nella componente iesima del vettore caratteristico).
Es., il vettore v=[1,1,0,1,0,0] e’ il vettore caratteristico del sottoinsieme
formato dal primo, secondo e quarto oggetto, rispettivamente di lunghezza 1,
2, e 4. Il pattern associato e’ quello in figura.
1
4
2
Esercizio: determinare il vettore
caratteristico degli altri pattern raffigurati
1
6
2
2
4
4
3
4
7
2
1
3
6
2
3
4
1
Fomulazione di set partitioning
gap garantito 2!!
Un pattern di impaccamento
è ammissibile
se la somma dei pesi degli oggetti
selezionati nel pattern è  Q.
Considera l’insieme R di tutti i pattern
ammissibili descritti da vettori binari vr
di n elementi 0/1, in cui vri=1 significa che
l’oggetto i è selezionato nel pattern r.
Formulazione PLI (2):
un’altra applicazione del modello di
SET PARTITIONING:
Min rR zr
rR zr vir = 1
zr  {0,1}
 i  I pi variabile duale
rR
Obiettivo:
risolvere in modo efficiente il rilassamento continuo di
PLI(2) (0zr1).
Problema:
Il numero di variabili (pattern di taglio o cutting pattern) e’
esponenziale.
La formulazione completa del modello PLI2 ha un numero
esponenziale di colonne
Procedimento iterativo
All’iterazione h si lavora su un sottoinsieme delle colonne
RhR e si risolve il rilassamento continuo del modello di
Set Partitioning con le sole colonne in Rh.
(MASTER PROBLEM)
Ogni vincolo iesimo ha associata una variabile duale pi.
Sia Ph=[pih] il vettore delle variabili duali ottime alla fine della
h-esima iterazione

per i pattern attivi in Rh, il costo ridotto delle variabili zr è non
negativo (cond ottimalita’ per il problema di minimo).
Il costo ridotto della variabile zr associata alla colonna vr è
1 - iI pih vir  0.
Se anche il costo ridotto delle variabili in RRh
(non esplicitate) fosse  0
allora la soluzione corrente con le sole colonne
in Rh sarebbe ottima per l’intero problema
definito su tutte le colonne in R.

Occorre individuare se esistono in RRh colonne
con costo ridotto <0…
…. senza dover enumerare esplicitamente tutti
gli elementi dell’insieme R!
Supponiamo di disporre di un oracolo che determina
la colonna di costo ridotto minimo
(in funzione del valore delle variabili duali correnti)
• Se la colonna di costo ridotto minimo ha costo ridotto negativo, non sono
soddisfatte le condizioni di ottimalita’ della programmazione lineare, e
tale colonna va aggiunta a Rh per migliorare la soluzione corrente.
Dopo aver aggiunto la nuova variabile si riottimizza, effettuando una
iterazione del simplesso primale.
• altrimenti (la colonna di costo ridotto minimo ha costo ridotto non
negativo, ) l’oracolo fornisce la dimostrazione di ottimalità della soluzione
corrente rispetto a tutto l’insieme R, anche se le variabili in RRh non
sono mai stati esplicitamente generate.
Infatti e’ come se tutte le colonne fossero presenti nel modello ma poste
a 0 (quindi fuori base) in quanto il loro costo ridotto risulta 0, essendo
per def maggiore del costo ridotto minimo il quale per hp e’ 0, e quindi
soddisfacendo le condizioni di ottimalita’ della programmazione lineare
per un problema di minimo sull’intero insieme delle colonne R..
Questo passo (l’oracolo) e’ formulabile come
un nuovo problema di ottimizzazione
noto come “subproblem” o “pricing” poiche’ determina
il valore del costo ridotto minimo e quindi anche un limite
inferiore al costo ridotto di tutte le colonne
• Funzione obiettivo:
minimizzare il costo ridotto
• Vincoli di ammissibilita’:
essere una colonna della matrice dei vincoli del modello di
Set Partitioning

essere il vettore caratteristico (a componenti 0/1)
di un pattern di packing ammissibile
Struttura del generatore di colonne:
(SUBPROBLEM)
Min (1-  iI pih vi ) = 1 + Min -  iI pih vi = 1 + Max  iI pih vi tale che
le variabili vi siano le componenti del vettore caratteristico di un pattern ammissibile 
valori 0/1 per cui valga  iI pi vi  Q
Fatti importanti:
nel subproblem, le vi sono le incognite del modello, mentre le pih sono i dati dell’istanza,
essendo note in quanto sono le variabili duali ottime calcolate alla iterazione h-esima del
master problem (comando getdual(nome-vincolo) di Mosel).
Per analogia si pensi al ruolo delle variabili x* - soluzione ottima del rilassamento
continuo di un’istanza del problema di zaino binario - nella procedura di determinazione
della cover + violata (problema di separazione).
Anche in tal caso le x* sono utilizzate nel sottoproblema di separazione come dei dati di
input del problema (valori noti) e non come variabili.
Le variabili del problema di separazione sono associate alla decisione “l’oggetto iesimo
appartiene alla cover cercata”?
Questa coincidenza non e’ casuale. Infatti il ruolo che gioca il problema di separazione
nel primale (il rilassamento continuo del problema di PLI in esame) e’ analogo al
problema di generazione di colonne del duale associato.
NB il problema di generazione di colonne e’ anche detto problema di “pricing” poiche’ “da un
prezzo” alle colonne in funzione del loro “costo ridotto”
Calcolare il pattern di costo ridotto minimo in R
equivale a risolvere un problema di zaino
Max  iI pih vi
 iI pi vi  Q
vi  {0,1}
iI
Sia v* la soluzione ottima:
• se il valore ottimo della funzione obiettivo è 1
( iI pihv*i 1)
allora la soluzione corrente è ottima (STOP),
• altrimenti si include il nuovo pattern in Rh
(Rh+1 = Rhv*) e si ITERA il procedimento sul nuovo
insieme di variabili Rh+1.
esercizi
1)
Implementare con Mosel la procedura appena descritta
(esercizio svolto in laboratorio, per il caso + generale in cui gli
oggetti hanno domanda intera  il sottoproblema di generazione
di colonne e’ uno zaino a variabili intere)
2) Descrivere in quale tipo di problema di ottimizzazione consista il
problema di generazione di colonne (pricing) per risolvere
attraverso la generazione di colonne il rilassamento continuo della
formulazione come problema di set partitioning del problema di
Airline Crew Scheduling
Generazione di colonne e
decomposizione di Dantzig-Wolfe
•
•
Dalla programmazione lineare ricordiamo che i punti di un politopo P corrispondenti alle
soluzioni di un sistema di vincoli lineari {x: Dxd} sono esprimibili come
somma della combinazione convessa dei vertici (xv V) di P
piu’ la combinazione conica dei raggi estremi (rr R) di P (quando P non e’ finito)
P = {x: Dx  d} = {x: x = vV xv lv + rR rr ar tale che vV lv =1, lv 0, ar 0 }
(cosidetta decomposizione di Dantig-Wolfe)
Si osservi che nella prima formulazione le incognite sono le x, mentre nella seconda le incognite
sono i coefficienti l e a, mentre i vettori che descrivono i vertici e i raggi del politopo sono
(teoricamente) noti.
Supponiamo di avere un problema del tipo Min {x: Ax b, Dx d, x Zn} e che
Min {x: Dx d, x Zn} sia “trattabile” (facilmente risolvibile, ad esempio quando il politopo {x: Dx
d,} ha vertici interi)
In tal caso puo’ convenire utilizzare una diversa formulazione del problema, in cui sostituisco
alla variabile x, nell’ambito dei vincoli Ax b, la somma della combinazione convessa dei vertici
+ combinazione conica dei raggi estremi del politopo {x: Dx d, x Zn}.
Il nuovo modello ha un numero molto elevato di colonne (sono i vertici e raggi di {x: Dx d, x
Zn}) per cui viene risolto attraverso la tecnica di generazione di colonne.
L’utilizzo congiunto di queste due tecniche e’ molto diffuso perche’ permette di sfruttare la
struttura matematica dei sottoproblemi.
Ne vediamo un esempio nel caso del cammino di costo minimo di durata limitata
Si cerca l’ottimo intero in P1P2: (punti rossi
in figura), X(P1P2)
caso specifico: il politopo P1 ha vertici non necessariamente interi,
mentre
Il politopo P2 ha vertici interi,quindi P2=CH(X(P2))
 Il politopo P3 dato dalla loro intersezione P1P2
ha alcuni vertici interi (quelli che sono anche vertici di P2) e altri no.
Vogliamo sfruttare il fatto che ciascuno dei vertici di P3, essendo
anche un punto interno di P2, e’ esprimibile come combinazione
convessa dei vertici di P2
x2
P2
P3
Indichiamo con X(P)
le soluzioni intere
appartenenti al politopo P,
e con CH(X(P)) la
combinazione convessa
dei punti interi in P
P1
x = l1x1 + l2x2
x1
In questo caso risolvere
il rilassamento continuo
di min cx: xP1P2
da lo stesso risultato
rispetto a risolvere
min cx: xP1CH(X(P2))
Si cerca l’ottimo intero in P1P2: (punti rossi
in figura), X(P1P2)
Caso generale: il politopo P1 ha vertici non necessariamente interi,
Il politopo P2 non ha vertici interi, ma e’ possibile calcolare la convex hull
dei suoi punti interi X che indichiamo come CH(X(P2))
 Il politopo P3 dato dalla loro intersezione P3 = P1CH(X(P2))
ha alcuni vertici interi (quelli che sono anche vertici di CH(X(P2)) ) e altri no.
Vogliamo sfruttare il fatto che ciascuno dei vertici di P3, essendo
anche un punto interno di P2, e’ esprimibile come combinazione
convessa dei vertici di CH(X(P2)), cioe’ di X(P2)
x2
CH(X(P2))
P3
In questo caso risolvere
il rilassamento continuo
di min cx: xP1P2
da un risultato peggiore
rispetto a risolvere
min cx: xP1CH(X(P2))
P1
x = l1x1 + l2x2
x1
P2
Cammino Minimo con Risorse Vincolate (RC-SPT)
Data una rete orientata e pesata G=(N,A) si cerca il cammino
di costo minimo dal nodo s al nodo t di durata inferiore a T
In assenza di altri vincoli (e di cicli di costo negativo sul grafo), il problema del
cammino di costo minimo (SPT) e’ facile e viene comunemente risolto da
algoritmi ad hoc, quali Dijkstra, Bellman-Ford, etc.
i
xij
j
In alternativa si formula un modello matematico come problema di flusso: una unita’
di flusso esce dal nodo s, origine del cammino, attraversa altri nodi del grafo e
giunge al nodo destinazione t.
Per formalizzare il percorso seguito dal flusso, si introducono delle variabili xij per
ogni arco (ij) in A, a valore compreso tra 0 e 1, e dei vincoli, detti vincoli di
bilancio, che impongono che
1) per ciascun nodo i,  da s e t., la differenza fra la somma dei flussi sugli archi
entranti nel nodo i e la somma dei flussi sugli archi uscenti da i sia nulla. Infatti
sia che il cammino passi o meno da i, in i esce la stessa quantita’ di flusso che
vi e’ entrata.
2) Per il nodo s, origine del cammino, la differenza e’ -1 (esce da s una unita’ e non
ne entra nessuna) mentre per t, destinazione del cammino, la differenza e’ +1
(entra una unita’ e ne esce 0).
Se un flusso xij attraversa l’arco (ij) paga un costo pari a xij cij.
Riflettiamo sulla forma che ha la matrice dei vincoli appena descritta (indicata con E):
Per ogni nodo esiste un vincolo di bilancio che impone un certo valore alla differenza fra flusso
entrante e flusso uscente dal nodo.
Quindi la matrice E ha una riga per ogni nodo del grafo e una colonna per ogni arco.
Poiche’ la variabile xij rappresenta il flusso presente sull’arco (i,j) allora la riga i-esima avra’
coefficiente -1 per ogni arco che esce dal nodo i e coefficiente +1 per ogni arco entrante
nel nodo i (0 per ogni arco che non e’ incidente nel nodo i).
i
7
s
5
2
La matrice di incidenza nodiarchi E
per il grafo orientato in figura
e’ data da
4
3
h
t
1
2
-1
(si) (su) (sj) (iv) (it) (ju) (uh) (uv) (hi) (ht) (vj) (vt)
6
s
i
j
u
h
v
t
-1
1
0
0
0
0
0
-1
0
0
1
0
0
0
-1
0
1
0
0
0
0
0 0
-1 -1
0 0
0 0
0 0
1 0
0 1
0
0
-1
1
0
0
0
0
0
0
-1
1
0
0
0
0
0
-1
0
1
0
0
1
0
0
-1
0
0
0
0
0
0
-1
0
1
0
0
1
0
0
-1
0
0
0
0
0
0
-1
1
v
6
j
8
1
u
In E, c’e’ una colonna
per ogni arco (i,j), con
2 sole componenti 0:
-1 nella riga i
+1 nella riga j
Il problema di cammino minimo, formulato come
problema di flusso, ha 12 variabili e 7 vincoli.
La funzione obiettivo Min CX e’ il prodotto scalare
tra il vettore delle variabili X e il vettore dei costi C
i
7
s
5
2
4
3
h
t
1
2
X = [x(si) x(su) x(sj) x(iv) x(it) x(ju) x(uh) x(uv) x(hi) x(ht) x(vj) x(vt)]
-1
C=[7
4
2
3
5
6
6
1
2
1 -1
8]
Per il problema di cammino minimo da s a t,
la formulazione estesa e’ la seguente:
-xsi - xsu -xsj = -1
-xiv - xit +xhi + xsi = 0
b e’ il vettore
xsj - xju +xvj = 0
dei bilanci ai nodi
xsu + xju -xuh - xuv = 0
xuh -xhi -xht = 0
xiv +xuv -xvj -xvt = 0
xit + xht +xvt = +1
e i vincoli di non negativita’ xij 0  (i,j) ∊ A
v
6
j
6
8
1
u
In forma compatta,
Min cx: Ex=b, x0
dove E e’ la matrice
di incidenza nodi/archi
La matrice E e’
TOTALMENTE
UNIMODULARE
Cammino Minimo con Risorse Vincolate (RC-SPT)
La matrice E dei soli vincoli di bilancio del flusso e’ unimodulare.
Tale proprieta’ rende il modello una formulazione ideale del problema,
cioe’ il rilassamento continuo (0 ≤ xij ≤ 1) ha soluzione ottima intera.
Infatti SPT e’ un problema “facile” (polinomiale).
Utilizzare il modello matematico per il problema del cammino di costo minimo (SPT) e’ utile
solo quanto sono presenti restrizioni aggiuntive (altrimenti si usano gli algortimi ad hoc).
Queste restrizioni sono facilmente formalizzabili nell’ambito del modello matematico,
aggiungendo ai vincoli di bilancio del flusso Ex=B, altri vincoli.
Al contrario, difficilmente gli algoritmi ad hoc sviluppati per SPT sono in grado di gestire al loro
interno vincoli addizionali.
L’esempio considerato in questo caso riguarda un attributo legato al tempo di percorrenza di
un arco, e un limite alla durata totale del cammino.
I cammini ammissibili hanno durata non superiore a T.
Il nuovo vincolo distrugge la unimodularita’ della matrice per cui il rilassamento continuo del
modello ha soluzione ottima frazionaria (vedremo perche’ con un semplice esempio).
Vediamo come cambia il modello precedente e studiamo un approccio di soluzione basato
sulla generazione di colonne e la DW decomposition per la soluzione del rilassamento
continuo.
Gli algoritmi ad hoc per SPT possono essere generalizzati per gestire questo specifico vincolo
(durata massima) ma la loro complessita’ diventa da polinomiale a pseudo-polinomiale e
la loro esecuzione tipicamente richiede molto tempo di calcolo.
Cammino Minimo con Risorse Vincolate (RCSPT)
modello per archi
Min Σ(ij)A cij xij :
obiettivo
1.
2.
3.
Σ(ij)FS(i) xij - Σ(hi)BS(i) xij =0  is,t
vincoli di bilancio Ex=b:
Σ(sj)FS(s) xsj = 1
garantiscono che gli archi selezionati
Σ(it)BS(t) xit = 1
siano un cammino da s a t
4.
Σ(ij)A tij xij ≤ T
5.
xij{0,1}
vincolo di durata max: distrugge la tot. unimodularita’
vincolo di interezza; adesso e’ necessario
La matrice dei soli vincoli di flusso 1-3 e’ unimodulare
MA
l’introduzione del vincolo 4 distrugge l’unimodularita’ della matrice dei vincoli 1-4, e rende
necessario il vincolo 5 per ottenere una formulazione intera.
La soluzione del rilassamento continuo del modello (1-5) e’ la combinazione convessa di
cammini di durata sia superiore che inferiore a T.
RC-SPT: arc based formulation
Min Σ(ij)A cij xij :
1.
2.
3.
Σ(ij)FS(i) xij = Σ(hi)BS(i) xij  i Perche’ la soluzione ottima del rilassamento
continuo e’ frazionaria?
Σ(sj)FS(s) xsj = 1
Il flusso uscente da s si distribuisce lungo diversi
Σ(it)BS(t) xit = 1
cammini, ciascuno di diverso costo e lunghezza.
4.
Σ(ij)A tij xij ≤ T
5.
xij{0,1}
Ciascun cammino concorre alla durata e al costo totale
In funzione della frazione di flusso che lo percorre
Esempio, ci siano solo due cammini disgiunti che collegano s a t,
il primo, p*, di costo c* e durata T*=T+D, quindi non ammissibile,
il secondo, p^, di costo c^ >c* e durata T^=T-d e quindi ammissibile e soluzione ottima intera
del problema.
La soluzione ottima del rilassamento continuo sara’ la combinazione convessa di p* e p^.
Si osservi che poiche’ sia p^ che p* rispettano i vincoli (1-3) anche ogni loro combinazione
convessa li soddisfa.
La generica soluzione sara’ (1-a)p^ + a p* con a  [0..1], e di costo (1-a)c^ + a c* < c^, e
tale che la durata (1-a)T^ + a T* soddisfi al vincolo (1-a)T^ + a T*  T.
Poiche’ al crescere di a diminuisce il costo della combinazione convessa, allora a sara’ il
massimo valore che consente alla combinazione convessa di rispettare il vincolo 4
sulla durata,
 a tale che (1-a)T^ + a T* = (1-a)(T-d) + a (T+D) = T + aD – (1-a)d = T

aD – (1-a)d = 0  a = d/(D+d).
Quindi lungo p^ passa flusso pari a (1-a) mentre lungo p* ne passa a.
Riformulazione per cammini (1)
Il poliedro associato ai vincoli (1-3) e (5) ha come matrice dei coefficienti la
matrice E e come vertici tutti i 9 camini p che vanno dal nodo s al nodo t.
Nel grafo in figura, con T=14, la soluzione ottima intera e’ il cammino
p3=(1,3,2,4,6) di costo 13 e durata 13.
Risolvendo con Xpress il rilassamento continuo, otteniamo la combinazione
convessa di p2 (l2=0.2) e p6 (l6=0.8) di costo 7 e durata 14
Esempio con T=14, s=1, t=6
Cammini (non) ammissibili, costo e durata:
11
2
1 10
4
12
1
17
2 3
10 3
10 1
57
3
6
22
5
12 3
cij tij
Sia P l’insieme di tali (9) cammini
p1
p2
p3
p4
p5
(1,3,2,4,5,6),
(1,3,2,5,6),
(1,3,2,4,6),
(1,2,4,5,6),
(1,3,4,5,6),
p6
p7
p8
p9
(1,2,5,6),
(1,2,4,6),
(1,3,5,6),
(1,3,4,6),
l2=0.2
l6=0.8
24,
15,
13,
14,
27,
9
10
13
14
13
5,
3,
24,
16,
15
18
8
17
model ModelName
uses "mmxprs";
Modello Mosel per la combinazione convessa di
costo minimo
e durata vincolata
dei cammini s-t
declarations
num_paths = 9
paths = 1..num_paths
costo : array(paths) of integer
tempo : array(paths) of integer
Tmax : integer
lambda : array (paths) of mpvar
Objective, durata:linctr
end-declarations
Tmax := 14
costo :: [24, 15, 13, 14, 27, 5, 3, 24, 16]
tempo :: [ 9, 10, 13, 14, 13, 15, 18, 8, 17]
forall (p in paths) lambda(p) <=1
sum(p in paths) lambda(p) = 1
durata := sum(p in paths) lambda(p)* tempo(p) <= Tmax
Objective := sum(p in paths) lambda(p)* costo(p)
minimize(Objective)
forall(p in paths| getsol(lambda(p))>0)
writeln("lambda(",p,")=", getsol(lambda(p)),
" costo=",costo(p)," tempo=",tempo(p))
write("costo ottimo continuo= ", getobjval,
"\ntempo ottimo continuo= ", getact(durata))
end-model
Output:
lambda(2)=0.2 costo=15 tempo=10
lambda(6)=0.8 costo=5 tempo=15
costo ottimo continuo= 7
tempo ottimo continuo= 14
x3
2
x1
4
x5
1
x9
x4
x2
x8
x6
3
x10
5
x7
6
Combinazione convessa dei cammini nello spazio c,t
T=14, s=1, t=6
Cammini,
costo durata
cammini
20
p7
18
p1
p2
p3
p4
p5
p6
tempo p7
p8
p9
p9
16
p6
tempo
14
p*
12
p3
ottimo intero
10
p1
p2
8
p5
p8
6
4
(1,3,2,4,5,6)
(1,3,2,5,6)
(1,3,2,4,6)
(1,2,4,5,6)
(1,3,4,5,6)
(1,2,5,6)
(1,2,4,6)
(1,3,5,6)
(1,3,4,6)
24
15
13
14
27
5
3
24
16
9
10
13
14
13
15
18
8
17
Min costo
2
0
0
5
10
15
costo
20
25
30
Si nota subito come p* =
(combinazione convessa
di p2 e p6 che interseca
il vincolo t=14) sia l’ottimo
Riformulazione per cammini (2)
Il generico cammino p e’ identificato dal vettore caratteristico yp nello spazio degli archi {0,1}n con
n=|A|=10 numero di archi, con ypij =1 se l’arco (i,j)p.
Ad esempio, il cammino p1 dato dalla sequenza di nodi (1,3,2,4,5,6) e’ composto dagli archi
(1,3),(3,2),(2,4),(4,5),(5,6).
Assegnando un ordine agli archi come quello in figura (arco1=(1,2), arco2=(1,3), arco3 = (2,4),
arco4=(2,5), arco5=(3,2), …., arco10=(5,6)),
il vettore caratteristico del cammino p1 e’y1 = [ 0, 1, 1, 0, 1, 0, 0, 1, 0, 1 ],
quello del cammino p2 e’ y2 = [ 0, 1, 0, 1, 1, 0, 0, 0, 0, 1 ],
e quello del cammino p6 e’ y6 = [ 1, 0, 0, 1, 0, 0, 0, 0, 0, 1 ].
…in modo analogo a come descriviamo i cutting pattern rispetto agli oggetti da tagliare!!....
x3
2
x1
4
x4
x9
x5
1
x8
x6
x2
3
x10
5
x7
6
p1
p2
p3
p4
p5
p6
p7
p8
p9
(1,3,2,4,5,6),
(1,3,2,5,6),
(1,3,2,4,6),
(1,2,4,5,6),
(1,3,4,5,6),
(1,2,5,6),
(1,2,4,6),
(1,3,5,6),
(1,3,4,6),
Riformulazione per cammini (3)
Il modello per cammini si ottiene dal precedente modello per archi, operando la
sostituzione xij =ΣpP lp yijp nel vincolo 4 e imponendo la combinazione convessa in P.
I vincoli 1-3 sono inglobati nei vettori y, poiche’ le y sono i vettori caratteristici dei vertici
del politopo descritto dai vincoli (1-3).
La soluzione del rilassamento continuo del modello Xpress per cammini, usa p2 con
moltiplicatore l2=0.2 e p6 con l6=0.8: manda 0.2 unita’ di flusso lungo p2 e 0.8 lungo p6.
Quindi sono coinvolti con flusso >0 gli archi di ciascuno dei due cammini, p2 e p6, mentre
tutti gli altri archi hanno tutti flusso nullo.
Nello specificio, x2 la variabile associata all’arco (1,3) ha flusso 0.2, mentre x10 la
variabile associata all’arco (5,6), ha valore 0.2+0.8=1, e x45=0 poiche’ (4,5)p2, p6.
x3
2
x1
4
x4
x9
x5
1
x8
x6
x2
3
6
x10
5
In generale l’arco (i,j) ha flusso
pari alla somma dei moltiplicatori
di tutti i cammini di cui (i,j) fa parte
(quelli nel cui vettore caratteristico
la componente associata
all’arco (i,j) ha valore 1
x7
Ogni soluzione del modello per cammini induce una e una sola
soluzione del modello per archi. Ex. Vale anche il viceversa?
p1
p2
p3
p4
p5
p6
p7
p8
p9
(1,3,2,4,5,6,
(1,3,2,5,6),
(1,3,2,4,6),
(1,2,4,5,6),
(1,3,4,5,6),
(1,2,5,6),
(1,2,4,6),
(1,3,5,6),
(1,3,4,6),
modello per cammini (1)
xij = ΣpP lp yijp
Min Σ(ij)A cij (ΣpP lp yijp) = ΣpP lp (Σ(ij)A cijyijp) = ΣpP
lp cp :
I vincoli 1-3 (che impongono che le x descrivano un cammino da s a t)
sono sostituiti da “pP” dove P e’ l’insieme dei cammini elementari da s a t
Σ(ij)A tij(ΣpP lpyijp) = ΣpP lp(Σ(ij)A tijyijp) = ΣpP lptp ≤ T
ΣpP lp yijp = xij
ΣpP lp = 1
lp  0
Le variabili di questo modello sono le xij e le lp
mentre le yijp sono coefficienti noti
xij{0,1}
Infatti, dato un cammino pP, e’ noto da quali archi e’ formato
(e’ scritto nel vettore caratteristico di ogni cammino p)
modello per cammini (2) rilassamento continuo
Min Σ(ij)A cij (ΣpP lp yijp) = ΣpP lp (Σ(ij)A cijyijp) = ΣpP lp cp :
1.
2.
3.
4.
5.
Σ(ij)A tij(ΣpP lpyijp) = ΣpP lp(Σ(ij)A tijyijp) = ΣpP lptp ≤ T
ΣpP lp = 1
ΣpP lp yijp = xij
p ≤ 1

0
≤
Σ
l
y
pP
p
ij
0 ≤ xij ≤ 1
lp  0

Nel modello rilassato, ricordando che i coefficienti yijp sono binari ({0,1}),
il vincolo di combinazione convessa
ΣpP lp = 1, e di segno lp  0
implica il vincolo 0 ≤ ΣpP lpyijp ≤ 1 che equivale ai vincoli [3,4].
Di conseguenza i vincoli [3,4] possono essere rimossi
Resta un problema di PL con un numero molto elevato di colonne
Infatti sono una per ogni cammino da s a t
modello per cammini (3):
Min ΣpP lp cp :
1. ΣpP lptp ≤ T
2. ΣpP lp = 1
3. lp  0
pP
(problema P)
t
a
Variabili duali
Come gestire la cardinalita’ esponenziale dell’insieme dei cammini P?
Procedimento iterativo, alla generica iterazione k :
Siano tk e ak i moltiplicatori duali ottimi associati ai vincoli 1 (t) e 2 (a) ottenuti
dalla soluzione del problema P rispetto all’ insieme corrente di cammini PkP.
Poiche’ tp= Σ(ij)A tijyijp e cp= Σ(ij)A cijyijp
il costo ridotto di un generico cammino pP all’iterazione k e’ dato da
ĉk,p = Σ(ij)A (cij - tk tij) yijp – ak.
Reinterpretando sul grafo G=(N,A),
certificare l’ottimalita’ della soluzione corrente rispetto a tutto P
equivale a dimostrare che il costo del cammino di costo minimo rispetto ai
coefficienti di costo ĉij k = cij - tk tij , e’  ak
Generazione di colonne:
calcolo di cammini minimi sul grafo G
secondo il vettore dei costi ĉ.
Ad ogni iterazione k il master problem risolve il problema di
programmazione lineare P nello spazio dei cammini correnti Pk e
ottiene i moltiplicatori duali tk e ak.
Il sottoproblema di generazione di colonne (pricing) risolve il problema di
cammino (elementare) di costo minimo da s a t nello spazio degli archi
rispetto ai costi modificati ĉk.
Sia p*k tale cammino e ĉ*k il suo costo.
Se ĉ*k ak allora la soluzione corrente e’ ottima rispetto a tutto P
altrimenti
Pk+1 = Pk  p*k (il cammino ottimo determinato viene aggiunto a Pk) e si
itera, risolvendo il master problem rispetto a Pk+1.
Punti critici del problema
di generazione di colonne per RC-SPT
Cicli di costo negativo:
Il grafo del nostro esempio e’ aciclico quindi ogni cammino s-t e’ elementare
(privo di cicli). Questo accade in tutti quei casi in cui il grafo modellizza delle
attivita’ definite nel tempo, e il cammino descrive sequenze di attivita’.
In caso contrario, il grafo puo’ contenere dei cicli. Se i costi ridotti degli archi (cij
- tk tij) diventano negativi, si possono formare cicli di costo negativo, in
presenza dei quali il problema di cammino minimo elementare diventa NPHard.
Ci sono diverse tecniche per ovviare a questo problema, ma resta un ostacolo
da tenere in considerazione nella scelta del metodo di soluzione se si lavora
con grafi non aciclici.
Velocita’ di convergenza
Per velocizzare la convergenza del metodo e’ utile aggiungere ad ogni
iterazione non solo la colonna (il cammino) di costo ridotto minimo ma tutti
quelli con costo ridotto <0.
Qualita’ del bound ottenuto risolvendo in Column
Generation il rilassamento continuo della
riformulazione secondo la DW decomposition.
•
Se P2 ha vertici interi il bound ottentuo e’ uguale a quello ottenibile dal rilassamento
continuo della formulazione originale.
Puo’ essere comunque vantaggioso adottare questa tecnica per sfruttare i vantaggi
computazionali dati dalla struttura di P2, quando i vincoli rispetto a cui si opera la DW
decomposition sono gli unici vincoli aggreganti di k sottoproblemi indipendenti.
Esempio: flusso capacitato indivisibile multicommodity di costo minimo.
•
Se P2 non ha vertici interi ma si sanno generare le colonne associate ai vertici della
Convex Hull (inviluppo convesso) dei punti interi di P2, allora il bound ottenuto e’
migliore di quello associato al rilassamento continuo della formulazione originale.
Vedremo nel seguito che questo bound eguaglia il valore del bound ottenibile
rilassando i vincoli “complicanti”.
•
In entrambi i casi la riformulazione spesso permette di gestire in modo efficiente la
generazione di colonne anche ai singoli nodi dell’albero di Branch and Bound (tecnica
di Branch and Price) velocizzando la soluzione del problema intero rispetto ad
applicare il Branch and Bound alla formulazione originale.
Esempio: flusso capacitato indivisibile
multicommodity di costo minimo
Data
• una rete G=(N,A) orientata, pesata (costo cij), e capacitata (capacita’ uij)
• K coppie {ok/dk } che sono origine/destinazione di una singola commodity di flusso
di domanda dk, k=1..K
Si vuole determinare come instradare tutto il flusso sulla rete in modo tale da:
• soddisfare la domanda, instradando ciascuna commodity lungo un singolo
cammino
• soddisfare i vincoli di capacita’ sugli archi
• minimizzare il costo totale
Questo problema rappresenta un caso esemplare rispetto alla possibilita’ data dalla
DW decomposition di calcolare il valore del rilassamento continuo del problema in
modo efficiente, sfruttando non solo la struttura del problema (problema di cammino
di costo minimo) ma anche la decomposizione in K sottoproblemi indipendenti una
volta rimosso il vincolo aggregante sulla capacita’ degli archi.
Inoltre, i cammini generati per ogni coppia o/d durante il processo risolutivo, sono
utilizzabili come input di un modello di PLI strutturato come il master problem che
fornisce una soluzione euristica al problema (non tutte le colonne sono presenti ma
solo quelle necessarie alla soluzione ottima del rilassamento continuo). Una tale
informazione non sarebbe desumibile dalla soluzione continua del modello per archi.
flusso capacitato indivisibile multicommodity di costo minimo
Sia E la matrice di incidenza nodi/archi del grafo G,
Per ogni commodity k
sia bk il vettore dei bilanci del flusso k-esimo,
quindi bkok= -1, bkdk= +1, bki= 0 per ogni nodo i  dk, ok.
Sia xkij la variabile booleana associata al transito del flusso di tipo k sull’arco (i,j)
e si indichi con xk il vettore dei flussi xkij di tutti gli archi (i,j)A
Sia xij la variabile intera associata al flusso totale che transita sull’arco (i,j)
Allora il problema puo’ essere modellizzato come il seguente problema di PLI:
Min ijA cij xij =  ijA cij (k dk xkij )
xij = k dk xkij  (i,j) A
! Il flusso di ogni k e’ indivisibile
xij  uij
 (i,j) A
! Vincoli di capacita’ su TUTTO il flusso in transito su (i,j)
E xk = bk
 k=1..K
!  k il flusso e’ bilanciato e la domanda e’ soddisfatta
xkij  {0,1}
 (i,j) A,  k=1..K
! se la commodity k usa o meno l’arco (i,j)
Gli unici vincoli aggreganti sono i vincoli di capacita’
Se non ci fossero, dovremmo solo risolvere K problemi di cammino di costo minimo da ciascun
nodo orgine ok alla propria destinazione dk (il flusso di ciascuna commodity e’ indivisibile quindi
segue lo stesso percorso)
flusso capacitato indivisibile multicommodity di costo minimo
Riformuliamo il problema secondo la DW decomposition applicata ai vincoli di flusso  k.
I vincoli di capacita’ sono gli unici vincoli in cui compare una sommatoria su k, rimossi i quali
il problema si decompone in k sottoprobemi indipendenti di cammino minimo (flusso indivisibile).
I vincoli di bilancio per ciascuna commodity si possono riformulare come la scelta di un
cammino dal nodo ok al nodo dk.
• sia P(k) l’insieme di tali cammini e sia yp{0,1}A con pP(k) il vettore caratteristico del
cammino p rispetto agli archi (i,j) in A (ypij =1 se l’arco (ij) sta nel cammino p, ypij =0 altrimenti)
• sia lp la variabile booleana associata alla scelta di tale cammino lungo cui instradare
tutto il flusso della commodity k.
• sia cp il costo del cammino p, allora la funzione obiettivo si puo’ riscrivere come:
Min  ijA cij xij =  ijA cij (k dk xkij ) = k dk  pP(k) lp cp
Si cerca un tale cammino  commodity k in modo che su nessun arco il flusso complessivo
violi la capacita’.
k dk pP(k) lp ypij  uij
 (i,j) A ! Vincoli di capacita’ sul flusso in transito su (i,j)
 pP(k) lp = 1
lp  {0,1}
 k=1..K
 p P(k),  k=1..K
! Un solo camino per ogni commodity
Poiche’ la cardinalita’ degli insiemi P(k) e’ molto elevata, si applica la column generation  k:
ad ogni passo si genera (se ha costo ridotto <0) un nuovo cammino per ogni commodity
risolvendo un problema di cammino minimo rispetto ai costi ridotti,
Utilizzando, per calcolare il costo ridotto degli archi, la variabili duali dei vincoli di capacita’.
Branch and Price:
integrazione del processo di generazione di colonne nel
Branch&Bound
•
Il procedimento di generazione di colonne serve per risolvere il rilassamento continuo di
modelli di PLI con un numero esponenziale di variabili.
•
Le colonne generate al nodo radice non necessariamente contengono tutte le colonne
presenti nella soluzione ottima intera. La soluzione ottima intera ottenuta rispetto ad esse
rappresenta una soluzione euristica subottimale.
•
Per determinare la soluzione ottima intera il procedimento di generazione di colonne va
integrato in un Branch and Bound.
•
Ad ogni nodo dell’albero di ricerca, la soluzione ottima del rilassamento continuo del
sottoproblema associato, puo’ richiede la generazione di ulteriori colonne.
•
La regola di branching utilizzata deve essere tale da essere facilmente integrabile nel
subproblem di pricing (si pensi alla citata regola del follow-on per air crew scheduling)
•
Tale procedimento prende il nome di BRANCH and PRICE.
La sua applicazione ha permesso negli ultimi 20 anni la soluzione di numerosissimi problemi di
scheduling, turnazione del personale, logistica e distribuzione per istanze reali non risolvibili
prima delle disponibilita’ di questa metodologia (turni dei piloti di AirFrance, scheduling dei
treni delle ferrovie olandesi, etc…). E’ indicativo quante volte questo procedimento sia
citato nelle applicazioni che hanno vinto il Franz Edelman Award nelle ultime edizioni
Rilassamento Lagrangeano
•
Esaminiamo tramite l’esempio del RC-SPT come operi il rilassamento Lagrangeano
e di quali proprieta’ goda.
•
Ricordiamo che la tecnica consiste nel rimuovere i vincoli “difficili” (la cui assenza
renderebbe “facile” il problema) dalla definizione della regione ammissibile,
trasferendo nella funzione obiettivo una componente che
– per i vincoli di disuguaglianza “premia” le soluzioni ammissibili (fornendo una valutazione
ottimistica del loro costo)
– per i vincoli di uguaglianza non altera il valore della funzione obiettivo delle soluzioni
ammissibili.
•
In entrambi i casi, ricordando la definizione di rilassamento, (allarga la regione
ammissibile o/e fornisce una valutazione ottimistica del costo delle soluzioni
ammissibili) si tratta di un rilassamento.
•
Operativamente, viene sommato alla funzione costo un termine dato dalla slack di
ogni vincolo moltiplicato per un peso detto moltiplicatore di Lagrange.
•
Quando il vincolo non e’ rispettato, la slack indica la violazione del vincolo stesso
•
Nel caso del RC-SPT il vincolo “complicante” e’ quello sulla durata massima del
cammino.
Rilassamento Lagrangeano del modello per archi del RC-SPT
P
Min Σ(ij)A cij xij :
Σ(ij)FS(i) xij - Σ(hi)BS(i) xij =0  is,t
Σ(sj)FS(s) xsj = 1
Si osservi che per ogni l 0 il problema P(l) rappresenta
un rilassamento di P. Infatti, per ogni x ammissibile,
Σ(it)BS(t) xit = 1
il termine Σ(ij)A tij xij - T e’ ≤ 0 e quindi la funzione obiettivo
Σ(ij)A tij xij ≤ T
fornisce una valutazione ottimistica del suo costo cx.
xij{0,1}
Sia x* la soluzione ottima di P, e sia x*(l) la
P(l)
soluzione ottima di P(l) per un dato l0.
Indichiamo con z* e F(l) il valore delle
Min Σ(ij)A cij xij + l (Σ(ij)A tij xij - T )
rispettive soluzioni ottime.
Σ(ij)FS(i) xij - Σ(hi)BS(i) xij =0  is,t Se x*(l)·t > T (quindi x*(l) non e’ ammissibile
Σ(sj)FS(s) xsj = 1
per P) allora, essendo x*(l) ottimo per P(l), vale
cx*(l) + x*(l)(t –T)  cx + l(tx –T) xP
Σ(it)BS(t) xit = 1
e quindi anche per il punto x* ottimo di P.
xij{0,1}
Da cio’ segue che
F(l) = cx*(l) + x*(l)(t –T)  cx* + l(tx* –T)  cx* = z*
Quindi  l0, F(l) e’ un lower bound di z*
Rilassamento Lagrangeano del modello per archi del RC-SPT
P(l)
Min Σ(ij)A cij xij + l (Σ(ij)A tij xij - T )
Σ(ij)FS(i) xij - Σ(hi)BS(i) xij =0  is,t
Σ(sj)FS(s) xsj = 1
Per ogni l 0 calcoliamo i costi ridotti secondo l: cij +l tij
Σ(it)BS(t) xit = 1
e risolviamo un problema di cammini minimo sul grafo G
secondo tali costi.
xij{0,1}
E’ come assegnare una penalizzazione a ciascun arco,
che si somma al costo vero e proprio, e che dipende dalla
durata del singolo arco secondo un peso l.
Vogliamo calcolare il peso “ottimo” l* in modo che la
combinazione convessa ottima dei cammini da s a t
che minimizza il costo ridotto sia la stessa che
risolverebbe la formulazione per cammini secondo DW
L’idea puo’ essere la seguente: partiamo da l1 (es, l1=0), risolviamo P(l1).
Se il cammino ottenuto supera la durata massima T, aumentiamo la penalita’ l
associata alla lunghezza degli archi nella formula del costo ridotto,
altrimenti proviamo a diminuirla.
Questo procedimento e’ alla base dell’ algoritmo del subgradiente.
Combinazione convessa dei cammini nello spazio c,t
ottimo continuo, intero e rilassamento lagrangeano
tempo
T=14, s=1, t=6
Cammini,
costo durata
p7
p9
T=14 p6
p*
p4
p3 ottimo intero
p2
p5
p1
p8
p1
p2
p3
p4
p5
p6
p7
p8
p9
(1,3,2,4,5,6)
(1,3,2,5,6)
(1,3,2,4,6)
(1,2,4,5,6)
(1,3,4,5,6)
(1,2,5,6)
(1,2,4,6)
(1,3,5,6)
(1,3,4,6)
24
15
13
14
27
5
3
24
16
9
10
13
14
13
15
18
8
17
Min costo
costo
Curva di livello della funzione obiettivo di P
Curva di livello del problema rilassato per
un certo l (funzione obiettivo min c+lt)
p*, combinazione convessa
di p2 e p6 che interseca
il vincolo t=14, e’ l’ottimo
continuo, e p3 quello intero
Quale Lambda? Una risposta la fornisce il rilassamento lagrangeano
Rilassamento lagrangeano e vincoli di uguaglianza
TH
Sia P = min {cx: Ax=b, xX} e sia P(l) = min {cx + l(b-Ax): xX}
Allora F(l) e’ un lower bound di z* per ogni lR.
Dim.
F(l) = cx*(l) + l(b-Ax*(l))  cx*(l) + l(b-Ax*(l)) = cx*
L’unica differenza tra il rilassamento di equazioni e quello di disequazioni e’ data dal vincolo di
segno del moltiplicatore associato al vincolo rilassato (in modo da penalizzare le x P).
Poiche’ F(l) fornisce un lower bound per P,
e’ ragionevole cercare quel valore di l che massimizza F(l), risolvendo il cosidetto
DUALE LAGRANGEANO:
Max F(l) tale che lR o lR+ a seconda dei casi (in caso di vincoli di disuguaglianza).
Nota che, se rilasso un set di n vincoli, lRn e’ un vettore con una componente per ciascun
vincolo rilassato. Il vettore l viene detto il vettore dei Moltiplicatori Lagrangeani.
Ogni moltiplicatore e’ vincolato in segno se tale vincolo e’ una disuguaglianza, con le stesse
regole che valgono per le variabili duali nella Programmazione Lineare, in modo da
penalizzare le x P .
Il duale lagrangeano
DL= Max F(l) tale che lR o lR+ a seconda dei vincoli rilassati (=, o  e
).
Interroghiamoci sulla forma della funzione F(l).
DL = max F(l) = maxl minx {cx + l(Ax-b): xX} =
al variare di l, che influenza solo la funzione obiettivo ma non i vincoli del problema di minimo, i
punti di ottimo x*(l) sono vertici della convex hull di X, indicata con CH(X).
Supponiamo di poter enumerare tali vertici nell’insieme {xv con vV, insieme dei vertici di CH(X)}
= maxl minvV {cxv + l(Axv-b): vV} =
Fissato il singolo xv, al variare di l l’espressione cxv + l(Axv-b) e’ lineare in l quindi identifica un
iperpiano in Rn (nello spazio dei l) 
DL consiste nel massimizzare nello spazio dei l una funzione concava poliedrale a tratti, tale che
per ogni l l’iperpiano attivo e’ quello associato alla soluzione ottima x*(l).
= maxl,wR { w: w  cxv + l(Axv-b)  vV } = maxl,wR { w: w - l(Axv-b)  cxv  vV } =
Questo e’ un problema di Programmazione Lineare nelle variabili w e l, con un numero elevato di
vincoli. Infatti c’e’ un vincolo per ogni vertice di CH(X).
Formuliamone il problema duale: c’e’ una variabile av  vincolo, vale a dire  vertice xv di CH(X)
= minaRV { vV av cxv : vV av =1, vV av (b-Axv) = 0, av 0  vV }
Oppure
minaRV { vV av cxv : vV av =1, vV av(b-Axv)0, av0  vV } se P = min {cx: Axb, xX}
che corrisponde alla riformulazione del problema P secondo DW rispetto ai vincoli Ax=b (Axb)
Il duale lagrangeano
DL= maxwR {w: w cxv + l(Axv-b)  vV }
DL consiste nel massimizzare nello spazio dei l una funzione concava
poliedrale (lineare a tratti), tale che per ogni l l’iperpiano attivo e’ quello
associato alla soluzione ottima x*(l)
Nei kinck points sono attivi piu’
iperpiani: il gradiente non e’ unico,
ma esistono infiniti subgradienti
F(l)
c x*(lj) + lj (Ax*(lh)-b)
Per ottimizzare questa funzione utilizziamo
l’algoritmo del subgradiente,
una generalizzazione per funzioni non differenziabili
del noto Algoritmo del Gradiente
l
lj
Costruiamo la funzione DL(l) per RC-SPT del nostro grafo
DL= maxwR {w: w cp + l(tp-14)  pP }
Ogni cammino identifica una retta che varia in funzione del valore di l, che possiamo
interpretare come il peso da attribuire alla durata unitaria degli archi
T=14
Cammini,
F(l)
c p tp
p11 (1,3,2,4,5,6) 24 9
p
15 10
p2
p2 (1,3,2,5,6)
3
(1,3,2,4,6)
p
p3
p4 (1,2,4,5,6)
p4
p5 (1,3,4,5,6)
(1,2,5,6)
p (1,2,4,6)
p6
p8 (1,3,5,6)
p9 (1,3,4,6)
p7
6
p
p5
7
1
2
l
3
6
p8
p9
13
14
27
5
3
24
16
13
14
13
15
18
8
17
p1
p2
p3
p4
p5
p6
p7
p8
p9
La funzione poliedrale F(l) e’ evidenziata in ROSSO
Algoritmo del sub gradiente
Supponiamo di risolvere DL= Max F(l) : lR (lR+ a seconda dei casi).
Nel caso di RC-SPT la funzione F(l) e’ minp in P cp+ l(tp-14)
Interpretiamo l come la penalizzazione dovuta alla durata.
Supponiamo di partire con l=0, e calcoliamo il cammino ottimo associato (p7).
Se p7 e’ ammissibile (t7-14)0, allora possiamo diminuire la penalita’, per permettere la
scelta di cammini + lunghi.
Al contrario, se p7 e’ non ammissibile, i.e. (t7-14)>0, il valore di l va incrementato per
scoraggiare l’utilizzo di cammini lunghi come e piu’ di p7.
Matematicamente, nell’espressione l(tp-14), la parte (tp-14) e’ il coefficiente di l in una
espressione lineare, quindi e’ il suo gradiente, ovverossia
la direzione locale di spostamento di massima crescita.
Poiche’ DL(l) e’ un problema di massimo, (tp-14) tale direzione e’ quella lungo cui
muoversi nell’aggiornare il valore di l (cerco la massima crescita)
Quindi lk+1 = lk + qk(tp-14)
In generale ci si muove lungo la direzione ( Ax*(l)-b )
Ora occorre definire una regola per la determinazione del passo qk
Algoritmo del sub gradiente:
determinazione del passo lungo la direzione di crescita (Ax*(l)-b )
Consideriamo il caso generale lk+1 = lk + qk(Ax*(l)-b ) = lk + qk(Axk -b)
dove per brevita’ xk indica la soluzione ottima per il vettore di moltiplicatori lk
Vogliamo determinare una regola per l’aggiornamento del passo qk
in modo tale da rispettare le seguenti
Condizioni necessarie per la convergenza del metodo: qk 0, k qk 
Se conoscessimo il valore ottimo di DL(l)=L*, vorremmo che il passo ci conducesse esattamente
in quel punto, vale a dire che
F(lk+1) = F( lk + qk(Axk-b)) = L*
Poiche’ lk+1 = lk + qk(Axk-b), e conoscendo l’espressione per F(l), vale che
F(lk+1) = cxk + [lk + qk(Axk-b)] (Axk-b) = L*
 un valore “ottimo” per lo step e’ qk = (L* - F(lk)) / (||Axk-b||2 )2
Quindi, la differenza tra il valore ottimo e il valore attuale (gap) della lagrangiana, normalizzato
per il quadrato della norma 2 del subgradiente nel punto corrente, e’ il valore ideale da dare
allo step. Non disponendo del dato L*, per valutare il gap possiamo approssimare L* col
valore di una soluzione ammissibile zH. eventualmente moltiplicato per un fattore correttivo fk
con 0 < fk < 2.
In conclusione, qk = fk (zH - F(lk)) / (||Axk-b||2 )2
Nota che se i vincoli rilasssati sono delle disuguaglianze per cui il moltiplicatore associato e’
vincolato in segno, la regola va generalizzata a
lk+1 = max { 0, lk + qk(Axk -b) }
Rilassamento Lagrangeano: condizioni di ottimalita’
• Il significato del subgradiente come direzione di crescita della
funzione lagrangeana, fornisce una condizione di ottimalita’ del
punto l: ||Ax*(l)-b||2 = 0
• Per il caso in cui i vincoli rilassati sono equazioni (Ax=b), se x*(l) e’
ammissibile per i vincoli rilassati e l e’ il punto di ottimo del duale
lagrangeano, allora x*(l) e’ ottimo.
• Se i vincoli rilassati in un problema di minimo sono del tipo Axb,
allora occorre anche che valgano le seguenti condizioni di
“complementarieta’”, (che nel caso della programmazione lineare
corrispondono agli scarti complementari, nell’ottimizazione non
lineare sono note come KarushKhunTucker conditions):
(Ax*(l)-b)·l = 0,
vale a dire che
- se il vincolo e’ rispettato da x come disuguaglianza stretta, allora il
moltiplicatore associato deve essere nullo,
- se un moltiplicatore li e’ >0, allora il vincolo associato, aix*(l)  bi
deve valere come uguaglianza.
Proprieta’ del rilassamento lagrangeano
Supponiamo di risolvere DL= Max F(l) : lR (lR+ a seconda dei casi).
Interroghiamoci sulla qualita’ del valore del rilassamento ottenuto rilassando in modo
Lagrangeano i vincoli Ax=b rispetto al valore ottenuto dal rilassamento continuo o a
quanto ottenibile tramite la DW decomposition sui vincoli xX.
Riscriviamo il problema P = min {cx: Ax=b, xX} come min {cx : Ax=b, Dx=d xZn} =
min {cx: x(XAXD)} dove XA e XD sono rispettivamente le soluzioni intere
contenute nei politopi associati ai vincoli Ax=b e Dx=d.
Abbiamo appena dimostrato che Max F(l) equivale a risolvere un rilassamento di P
ottenuto applicando la decomposizione di Dantzig-Wolfe sui vincoli xX, e quindi a
risolvere min {cx : Ax=b, xCH(XD)}
Quando risolvere min {cx : xXD} equivale a risolvere min {cx : xCH(X)} allora il
rilassamento lagrangeano rispetto ai vincoli Ax=b equivale a risolvere il
rilassamento continuo.
Non ci sono vantaggi dal punto di vista del bound ottenuto ma ci sono vantaggi
computazionali. Inoltre, la soluzione rilassata non ammissibile puo’ essere
facilmente “aggiustata” per rispettare i vincoli rilassati, dando luogo a euristiche
lagrangeane, che ad ogni iterazione calcolano sia il bound dato dal valore corrente
di F(l) che una soluzione primale ammissibile.
Infine, rispetto ad applicare la DW decomposition, risolvere il duale Lagrangeano e’
computazionalmente piu’ efficiente, ma per avere la soluzione primale continua
esatta occorrono algoritmi sofisticati per l’ottimizazione non lineare.
Il politopo del rilassamento continuo
{x: Axb  Dxd}
Le soluzioni intere del problema
{x: Axb, Dxd, xZ+}
Axb
Axb
4
Riferimento
Bibliografico:
Ahuja et al.
Network Flows
Cap 16
4
3
2
3
2
1
1
1
Axb
4
2
3
1
4 Dxd
La Convex Hull delle soluzioni
intere di {x: Dxd, xZ+}
Axb
4
2
3
4
L’insieme dei punti
{x: Axb  CH(Dxd, xZ+)}
3
3
2
2
1
1
1
2
3
4
1
2
3
4
Focus sui decisori
•
Fino ad ora abbiamo visto problemi in cui era presente un unico decisore con una unica
funzione obiettivo.
La realta’ e’ spesso ben piu’ complessa
•
Quando sono presenti diversi criteri di valutazione della soluzione, questi (tranne uno)
sono trattati come vincoli a cui si pone una soglia di accettazione: es, se un cammino
viene valutato sia secondo il costo che secondo la durata, si stabilisce una gerarchia fra
criteri (es il costo e’ + importante dela durata) e si fornisce una soglia di accettazione per
il valore del criterio meno importante (es. La durata del camino non deve superare Tmax)
e si ottimizza secodno il secondo criterio.
Questo e’ il caso piu’ comune ma non contempla numerose altre possibilita’
•
Cosa accade se 1) il decisore non ha una gerarchia fra i due criteri oppure 2) se esistono
piu’ decisori, ciascuno di pari importanza, quindi tutti allo stesso livello gerarchico, e
ciascuno vuole ottimizzare un proprio criterio?
 L’ottimizzazione multiobiettivo fornisce gli strumenti risolutivi per questi problemi.
•
Cosa accade se i decisori sono su scale gerarchiche diverse?
 L’ottimizzazione bi (o multi) livello bi(multi)-livello fornisce gli strumenti risolutivi per questi
problemi
Ottimizzaione multiobiettivo
•
Quando i diversi criteri non sono fra loro correlati, l’ottimo secondo ciascun obiettivo
e’ diverso dall’ottimo secondo gli altri obiettivi.
Ad esempio, se i collegamenti + veloci costano di piu’, il cammino di costo minimo
non sara’ il cammino piu’ veloce e viceversa.
In quei casi non esiste un singolo ottimo globale, cioe’ una singola soluzione che sia
ottima per tutti i criteri.
•
Si passa al concetto di OTTIMO PARETIANO
•
Informalmente
– un punto si dice Ottimo Paretiano se non e’ possibile migliorare il valore di un criterio senza
peggiorare il valore di un altro.
– Un punto si dice dominato da un altro punto se il vaore di tutti I criteri nel secondo punto e’
non peggiore rispetto al primo e migliore per almeno un criterio.
– Gli ottimi paretiani sono punti non dominati
– L’insieme dei punti non dominati forma la frontiera di Pareto. Esssa viene anche detta
“trade off surface” poiche’ rappresenta I migliori punti di compromesso tra gli attributi.
– Risolvere un problema multiobiettivo consiste nel calcolare la sua frontiera di Pareto
– Generalmente questa e’ molto estesa, per cui ci si limita a voler campionare tale frontiera
in un insieme di punti rappresentativi della frontiera stessa, non concentrati in una parte
ristretta, ma il piu’ possibile distanti fra loro.
– Si rivelano utili metodi risolutivi di tipo iterativo che esplorano la regione ammissibile
contemporaneamente in piu’ punti (metodi population based, quali gli algoritimi genetici),
all’interno dei quali si possa codificare un criterio di “diversita’ alle soluzioni gia viste.
Ottimalita’ e Dominanza
(caso di minimo di f1 e minimo di f2)
• Un vettore soluzione x1 domina x2 se
–
e’ strettamente migliore di per almeno un
obiettivo
– x1 e’ almeno tanto buono quanto x2 in tutti gli obiettivi
x1
x2
P = min f(x):
g(x) 0
xF
h(x)=0
x Rn, f(x) Rk
}
f2(x)
F
x1
x2
f1(x)
f2(x)
• Ottimalita’ Locale
–
F
e’ localmente ottimo se esiste un r>0 tale per cui
nessun punto x  FB(x1,r) (la palla di raggio r e
centro x1) tale che x domina x1. Ad esempio x2 non
domnina x1 poiche’ e’ peggio secondo f2.Il punto x3
domina x1 ma non appartiene all’intorno B(x1,r).
x1
x1
x2
x3
B(x1,r)
f1(x)
f2(x)
• Ottimalita’ Globale
F
– x e’ globalmente ottimo se non esiste alcun punto
nella regione ammissibile F che domini x.
x
f1(x)
• Cono negativo indotto da x1 in Rk
(spazio degli obiettivi)
P = min f(x):
g(x) 0
h(x)=0
x Rn, f(x) Rk
C- = {x: f(x) 0}
• Teorema del contatto
x1 e’ globalmente Pareto ottimo se
l’intersezione tra il cono negativo puntato
in x1 e la regione F coincide con x1.
(C- + x1)  F = {x1}
Se centriamo gli assi nella
soluzione A, allora A e’ dominata
da ogni soluzione nell’area
di dominanza,
e’ preferita rispetto a ogni soluzione
nell’area di preferenza, ed e’
incomparabile rispetto
a ogni soluzione
nelle aree di indifferenza.
f2(x)
x1
C-+x1
f1(x)
C-
f2(x)
Area di
Area di
indifferenza preferenza
CArea di
dominanza
A
F
Area di
indifferenza
f1(x)
Esempio
Consideriamo 5 punti soluzioni ammissibili del problema
max f1(x), min f2(x) con i valori di funzione obiettivo f1 e f2
in tabella
punto
f1
f2
A
8
5
B
9
2
C
12
1
D
11
2
E
16
2
f2(x)
A
5
4
3
B
D
E
2
Per ogni coppia di soluzioni X e Y
si puo’ riassumere la loro relazione
con una tupla di k=2 simboli,
+, -, =,
a seconda che X sia meglio, peggio
o uguale a Y per quell’ obiettivo
1
C
4
8
12
16
f1(x)
f2(x)
A
5
Max f1
Min f2
cono di
dominanza per A
4
3
B
D
E
2
f1
f2
A
8
5
B
9
2
C
12
1
D
11
2
E
16
2
C
1
4
A
A
8
12
B
C
D
E
–, –
–, –
–, –
–, –
–, –
–, =
–, =
+, +
–, +
B
+, +
C
+, +
+, +
D
+, +
+, =
–, –
E
+, +
+, =
+, –
–, =
+, =
La matrice e’ “anti-simmetrica”
Es. AB=(+,+)  BA=(–, –)
16
f1(x)
ESEMPIO…B e’ meglio di A
sia secondo f1 (massimizzo f1)
poiche’ f1(A) < f1(B)
( il 1 elemento della coppia AB e’ – )
sia secondo f2 (minimizzo f2)
poiche’ f2(A) > f2(B)
( il 2 elemento della coppia AB e’ – )
B cono di
dominanza di A

B domina A
f2(x)
A
5
r4
4
3
r3
B
r2
D
Max f1
Min f2
E r1
2
C
1
4
A
A
8
12
B
C
D
E
–, –
–, –
–, –
–, –
–, –
–, =
–, =
+, +
–, +
B
+, +
C
+, +
+, +
D
+, +
+, =
–, –
E
+, +
+, =
+, –
–, =
16
f1(x)
A domina ,
B domina {A}
C domina {A,B,D}
D domina {A,B}
E domina {A,B,D}
e’ dominato da {B,C,D,E}
e’ dominato da {D,C,E}
e’ dominato da 
e’ dominato da {C}
e’ dominato da 
+, =
Si dicono punti di rango 1 quelli non dominati da nessuno: punti C e E.
Sono di rango 2 quelli dominati solo da punti di rango 1: punto D.
Hanno rango 3 quelli dominati solo da punti di rango 1 o 2: punto B.
Il punto A e’ di rango 4.
Ottimizazione BiLivello
•
Quando esiste una relazione di gerarchia fra due (classi di) decisori,
l’ottimizzazione multiobiettivo non rappresenta correttamente il problema.
•
Infatti esiste un doppio livello decisionale, secondo cui il leader (decisore di primo
livello) impone “le regole”, e il follower (decisore di secondo livello) cerca di
ottimizzare la propria funzione obiettivo – diversa da quella del leaderrispettando le regole imposte dal leader.
•
Quindi ogni decisione del leader definisce la regione ammissibile del follower, che
deve rispettare le regole imposte dal leader. In questo insieme il follower ottimizza
la sua funzione obiettivo.
•
Quale e’ la regione ammissibile su cui il leader deve ottimizzare la propria
funzione obiettivo? A priori il Leader non sa quanto vale la sua scelta, poiche’
questo dipende dalla reazione del follower. Infatti il leader ha come uniche
possibilita’ le “risposte” del follower a ciascuna delle sue scelte. Si denota tale
insieme con eil “reaction set”.
•
Si osservi che la soluzione ottima (generalmente e’ unica a meno di soluzione
equivalenti) non necessariamente appartiene alla curva di Pareto del problema
che si otterrebbe mettendo tutti i decisori a pari livello, proprio perche’ in quel caso
viene a mancare una corretta rappresentazione della gerachia di decisione.
Uno strumento di prevenzione per progettare
in funzione del potenziale pericolo
•
Ci sono diversi punti in comune tra l’ottimizzazione bilivello e le applicazioni economiche
della teoria dei giochi, in particolare lo “Stackelber game” e l’equilibrio di Nash, da cui
provengono i concetti di leader e follower.
•
Quando il follower ha come obiettivo contrastare il leader (provocargli il maggior danno
possibile) l’ottimizzazione bilivello fornisce gli strumenti matematici per elaborare la migliore
strategia difensiva in termini di prevenzione, ovverossia individuare la scelta che minimizza il
peggior danno possibile, cioe’ quella che cerchera’ di porre in atto il follower. Questa filosofia
va sotto il nome di tecnica del minimax regret concept in decision analisys.
•
Dopo 9/11, negli USA sono state investire moltoe risorse per la protezione delle infrastruture
strategiche. Per quanto ingenti in ciascun progetto tali risorse sono limitate. Da cio’ sono nati
numerosi problemi reali legati al miglior utilizzo di tali risorse per “minimizzare il danno” in
caso di attacco.
•
La programmazione bilivello fornisce i modelli teorici e gli strumenti risolutivi adatti a questo
contesto.
•
In questa lezione ci limitiamo ad illustrare alcuni esempi di problemi di questo tipo e
vederemo un approccio risolutivo per una classe di problemi con una struttura matematica
specifica, basata sulle proprieta’ della programmazione lienare.
Per altre classi di problemi, occorre utilizzare strumenti tipici della Non Linear Optimization
(NLO), quali le condizioni di ottimo di Karush Khun Tucker (KKT).
•
Ottimizzazione bilivello e scelte economiche
I problemi bilivello trovano numerose applicazioni in campo economico, dove spesso
ci sono decisori con interessi contrapposti che appartengono a livelli gerarchici
diversi.
Una specifica area di utilizzo di questi modelli riguarda i problemi di tariffazione.
In mercati concorrenziali quali la telefonia e (soprattutto in altri paesi) l’energia, i
gestori di un servizio devono decidere a quale tariffa vendere i diversi prodotti
che offrono, in modo tale da ottimizzare il proprio profitto. E’ necessario tenere
conto del fatto che tale profitto non e’ costante poiche’ dipende dalla domanda.
A seconda del prezzo di ciascun servizio (decisione del leader) gli utenti (i follower)
modificano la loro domanda sia variando la quantita’ del singolo servizio che
passando da un servizio ad un altro con piano tariffario diverso, fino ad
abbandonare del tutto il gestore, rivolgendosi ai concorrenti.
Ci sono forti legami tra i concetti sviluppati nell’ambito dell’ottimizzazione bilivello e
le applicazioni della teoria dei giochi in campo economico, quali il gioco di Von
Stackelberg e l’equilibrio di Nash, da cui vengono mutuati i concetti di leader e di
follower, e di equilibrio.
Un punto di Nash-equilibrium (NE) e’ una soluzione rispetto alla quale nessuno dei
due giocatori ha interesse a modificare le proprie scelte. Questo e’ un punto di
ottimo del gioco di Stackelberg. Pero’ no ntutti I problemi godono della proprieta’
di esistenza di un punto di NE. Vediamone un esempio.
Localizzazione di servizi: maximum coverage
•
Nell’ambito dei problemi di “localizzazione dei servizi” esistono numerose varianti.
•
Una di queste consiste nel massimizzare il numero dei clienti serviti da p servizi che
devono essere localizzati in altrettante localita’ da scegliere nell’insieme dei clienti stessi,
con la regola che un servizio localizzato nel nodo j serve tutti i clienti posti entro una
distanza d da j.
•
Nella versione pesata ogni cliente ha un peso wi, e si vuole massimizzare la somma dei
pesi dei clienti che sono entro una distanza d dal piu’ vicino centro di servizio, attenuata
della distanza dalla facility. Nella versione non pesata si massimizza il numero dei clienti
serviti.
•
Entrambe le versioni di questo problema sono note come problemi di maximum converage,
evengono utilizzate per localizzare servizi pubblici.
•
Nella variante bilivello, il leader localizza i servizi ma ha un avversario che puo’ distruggere
un singolo link tra cliente e facility. L’avversario, che incarna il follower, ha come obiettivo
quello di procurare il maggior danno possibile al leader. Supponiamo che il follower
distrugga il link (i,j) tra il nodo cliente i e la facility j. Il danno provocato puo’ essere dovuto
al mancato servizio del cliente i collegato dal link distrutto all’unica facility j entro la distanza
d, se tutte la altre facilities sono state localizzate a distanza >d, oppure e’ dovuto alla
riduzione della qualita’ del servizio (nella versione pesata del problema) se il cliente e’ a
distanza <d da un’altra facility piu’ lontana.
Localizzazione di servizi: maximum coverage
•
•
•
Possiamo costruire delle istanze del porblema di
maximum coverage per cui non esiste un equilibrio di
Nash, cioe’ nelle quali a turno ogni giocatore (il
leader e il follower) hanno interesse a cambiare la
propria scelta.
Consideriamo per esempio 4 nodi posti ai vertici di
un quadrato di lato 1 e sia d=1.
Se il leader puo’ localizzare solo p=1 facility allora
sceglie un nodo a caso, es il nodo 1.
In questo modo serve se stesso e gli altri due nodi
adiacenti. Il follower distrugge uno dei due lati
adiacente al nodo 1, riducendo a 2 il numero dei nodi
serviti.
Allora il leader cambia il nodo in cui localizzare il
servizio e lo sposta sul nodo 3.
A sua volta il follower distruggera’ uno dei link
incidenti in 3, e cosi’ via.
Se invece il leader potesse localizzare p=2 facilities e
il follower distruggere, come prima, un unico link,
ponendo le faciliies su nodi non adiacenti si h un
punto di Equilibrio di Nash, che e’ anche l’ottimo del
problema
2
1
l =1
√2
4
3
1
2
4
3
Metodi di soluzione per i casi in cui il problema di ottimizzazione
affrontato dal follower e’ un problema di porgrammazione lineare
(oppure e’ un problema facile di PLI)
•
La parte complicata del modello matematico e’ scrivere il fatto che la funzione
obiettivo dipende sia dal valore delle variabili decisionali del leader (x) che da
quelle del follower (y), ma queste ultime a loro volta dipendono dalle scelte del
leader (x).
•
Indichiamo come y(x*) la soluzione y che ottimizza la funzione obiettivo del
follower quando x* e’ la decisione del leader, e indichiamo con xX e yY(x) i
rispettivi vincoli di ammissibilita’ delle due variabili.
Si osservi che poiche’ la regione ammissibile del follower e’ generalmente
determinata dalle decisioni y del leader, indichiamo questa relazione con Y(x).
•
Max cx + fy:
xX
y = argmin {dy : y Y(x) }
Esempio numerico
3,6
6
Min F(xy) = (x - 4y)
tale che:
x0
y = argmin f(y) = y:
x+y  3
2x – y  0
2x + y  12
3x – 2y  4
y0
y
5
Il punto di ottimo
(4,4) non e’
Pareto Ottimo!!
Regione
ammissibile
di y per x=3
4,4
4
gradienti
di F(x,y)
e di f(y)
3
Solo una parte della frontiera
2
del politopo costituisce
la regione ammissibile di x, in cui.
l’unico punto di ottimo e’ (4,4)
1
Se F(x,y) fosse stata da
massimizzare, avremmo avuto
due ottimi locali, (2,1) e (4,4).
2,1
Reaction set
1,2
1
2
Cono polare di -F e –f
Contiene i punti che
dominano (4,4)
secondo Pareto
3
4
x
5
Metodi risolutivi:
riformulare max/min f(y) tramite le condizioni di ottimalita’ di f(y):
un esempio, Guardie a Ladri
•
Una banda di contrabbandieri porta sigarette da s a t lungo il percorso meno
pericoloso.
•
La polizia cerca di contrastare il fenomeno pattugliando i singoli tratti da una
localita’ all’altra.
•
Se un tratto e’ pattugliato, i contrabbandieri non possono percorrerlo, altrimenti
verrebbero intercettati.
•
La polizia pero’ ha solo 3 pattuglie.
•
Obiettivo dei contrabbandieri: data la dislocazione delle pattuglie, seguire il
percorso + sicuro tra quelli percorribili.
•
Obiettivo della polizia: disporre le 3 pattuglie in modo da contrastare il piu’
possibile l’attivita’ dei contrabbandieri  rendere meno sicuro possibile il cammino
piu’ sicuro
Possiamo dare un modello matematico al
problema di guardie (leader) e ladri (follower)?
Sia ca il rischio di essere catturati lungo l’arco a.
I contrabbandieri cercano il cammino da s a t meno rischioso, risolvendo
Min cx : Ex=b, x0
Sia za la variabile booleana che descrive se l’arco a e’ sorvegliato.
La polizia dispone di 3 pattuglie. a za  3
Se l’arco a e’ sorvegliato (za=1) allora non e’ percorribile (xa=0), il che viene
espresso dal vincolo xa  1-za
Lo scopo della polizia e’ massimizzare il rischio del cammino scelto dai
contrabbandieri.
Sia w un lower bound alla sua pericolosita’
w  cx.
Allora l’obiettivo della polizia e’ max w
Un modello matematico bilivello
al problema di guardie e ladri
Max w:
w  cx
xa  1-za a A
a za  3
za {0,1} a A
OUTER
Problem,
Leader
(polizia)
x tale che x = argmin cx:
Ex=b
xa {0,1} a A
INNER
Problem,
Follower
(ladri)
dove
E e’ la matrice di incidenza nodi archi del grafo G=(N,A)
b e’ il vettore dei bilanci ai nodi che descrive il flusso di una unita’ dal nodo s al nodo
t.
Possiamo ricondurci a un problema con una sola funzione obiettivo?
Sfruttiamo le caratteristiche del inner problem (quello del follower)
Il duale del problema di cammino minimo
Il problema di camino minimo, min cx: Ex=b, xa{0,1} e’ equivalente al suo rilassamento
continuo min cx: Ex=b, xa0, per la totale unimodularita’ della matrice di incidenza E.
Quindi ammette un problema duale (detto problema di POTENZIALE) che ha lo stesso valore
di funzione obiettivo nel punto di ottimo (Th forte di dualita’)
Le variabili duali sono definite sui nodi, pi i N, non vincolate in segno.
Ricordiamo che il vettore dei bilanci ai nodi b e’ -1 per i=s, +1 per i=t, 0 else.
La funzione obiettivo del duale e’ quindi max pt - ps
Poiche’ la matrice E non ha rango pieno (ha rango n-1), si puo’ settare a un valore qualsiasi
una delle variabili duali ad esempio ps =0, in modo da avere una unica soluzione ottima
duale.
La funzione obiettivo si riduce a max pt
Nel primale c’e’ una variabile (colonna di E) per ogni arco (ij), quindi nel duale c’e’ un vincolo
per ogni arco. Ogni colonna (ij) di E ha -1 nella riga i, +1 nella riga j e il resto 0.
Max pt :
pj - pi  cij  a A, con a=(i,j)
Dalla dualita’ sappiamo che, se x e’ una soluzione ammissibile per il problema di cammino
minimo e p lo e’ per il problema di potenziale e inoltre cx=pt,
allora x e p sono OTTIME
Un modello matematico singolo livello
per il problema di guardie e ladri
Max w:
w  cx
xa  1-za
a za  3
za {0,1}
a A
a A
Ex=b
pj - pi  cij
cx = pt
 a A, con a=(i,j)
xa 0
a A
era
x = argmin cx:
Ex=b
Xa 0
Abbiamo introdotto le variabili e i vincoli del problema duale dell’inner problem, e abbiamo
sostituito alla richiesta che x fosse soluzione ottima del primale il vincolo di uguaglianza tra la
funzione obiettivo del primale e quella del duale.
Questa operazione e’ corretta in quanto l’inner problem e’ un problema di programmazione
lineare. Per altri modelli esistono tecniche simili ma + complesse basate sulle condizione di
ottimalita’ non lineari di KKT
Matrici (totalmente) Unimodulari
Def: A  Zmxn si dice Unimodulare 
ogni sua sottomatrice quadrata di rango
massimo ha determinante +1 o -1
Def: A  Zmxn si dice Totalmente Unimodulare
(TUM)  ogni sua sottomatrice quadrata ha
determinante +1 o -1
TH:
se A e b sono interi e A e’ unimodulare allora il
sistema Ax=b x0 ha soluzioni intere
TH:
se A e b sono interi e A e’ totalmente
unimodulare allora i sistemi
Ax  b x0 e Ax  b x0
hanno soluzioni intere
Molti problemi su reti di flusso godono di
questa proprieta’
Una matrice A con elementi 0, ±1 è TUM se
valgono le seguenti condizioni: (cond suff)
a) ogni colonna contiene al più due elementi non
nulli;
b) gli indici {1,…,m} delle righe di A possono
essere partizionati in due insiemi tali che ogni
colonna con i due elementi non nulli di segno:
- uguale, allora gli indici delle righe
corrispondenti appartengono a insiemi diversi;
- opposto, allora gli indici delle righe
corrispondenti appartengono allo stesso
insieme.
esercizio
Ricordiamo:
sia G=(N,A) un grafo orientato di n nodi N={1..n} e
m archi A={(i,j)  arco da i a j},
allora
la matrice E ha una riga iN e una colonna (ij)A.
La colonna (i,j) ha valori:
0 in ogni componente i,j
-1 nella componente i,
+1 nella componente j.
•
Verificare se la matrice E di incidenza nodi-archi di un grafo orientato
e’ TUM.
Per grafi non orientati E ha solo componenti +1 e 0.
•
Cosa accade per grafi non orientati?
esempio
sia G=(N,A) il grafo orientato in figura, con n=6, m=8.
Costruiamo la matrice E, di 6 righe e 8 colonne
La colonna (i,j) ha valori:
0 in ogni componente i,j
-1 nella componente i,
+1 nella componente j.
(12) (13) (24) (35) (43) (46) (52) (56)
2
1
2
3
4
5
6
-1
+1
4
-1
-1
+1
+1
-1
+1
+1
-1
-1
+1
+1
6
1
-1 -1
+1
3
5
esempio
Ma se G=(N,A) non fosse orientato ?
Costruiamo la matrice E, di 6 righe e 8 colonne
La colonna (i,j) ha valori:
0 in ogni componente i,j
+1 nella componente i,
+1 nella componente j.
2
4
(12) (13) (24) (35) (43) (46) (52) (56)
1
2
3
4
5
6
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1 +1
+1
6
1
3
5
La matrice E e’ ancora totalmente unimodulare.
E’ vero per tutti i grafi non orientati o esistono casi specifici in cui E non e’ TUM?
Di che proprieta’ gode il grafo G?
Proviamo a fare una visita in ampiezza da 1 colorando in modo alterno i nodi visitiati
a livello pari e a livello dispari
esempio
visita in ampiezza da 1 colorando in modo alterno i nodi visitiati a livello pari e a livello dispari
La procedura da lo stesso risultato a partire da ogni nodo!!
2
4
(12) (13) (24) (35) (43) (46) (52) (56)
1
2
3
4
5
6
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1 +1
+1
6
1
3
5
Grafi con questa propieta’ si dicono BIPARTITI
Una definizione analoga e’ che contengono solo cicli con numero pari di archi
Accoppiamento di massima cardinalita’
su grafo bipartito
G=(N,A) e’ bipartito

N=N1 N2
A  N1N2
Gli archi descrivono le
possibili coppie (matching)
Quante coppie al piu’?
Max ijA xij:
jN2 xij1 iN1
iN1 xij1 jN2,
xij{0,1}
esercizio
• Scrivere la matrice di incidenza del grafo in
figura, numerando da 1 a n1 i nodi in N1 e da
(n1+1) a (n1+ n2) i nodi in N2.
1
2
6
7
8
3
9
4
10
5
11
12
La matrice E dei vincoli e’ TUM
• Dimostrazione: le colonne soddisfano alla condizione (a) (due soli
elementi  0) e alla condizione (b) se partizioniamo gli indici delle
righe secondo la partizione dei nodi corrispondenti nei due
sottoinisiemi di nodi N1 e N2.
• Posso risolvere il problema con la programmazione lineare 
l’accoppiamento di massima cardinalita’ su grafi bipartiti e’ un
problema FACILE
• Esistono algoritmi polinomiali appositi (+ efficienti del simplesso)
• ? Quanto incide il fatto che il grafo sia bipartito? ....
Th Konig
Se G e’ bipartito allora 
Il matching di massima cardinalita’ ha lo
stesso valore del vertex cover minimo
(total dual integrality)
Dim:
sapendo che la matrice dei vincoli e’ TUM, scrivere il
duale del problema di matching, verificare se la nuova
matrice e’ TUM, e applicare il teorema forte di dualita’
Oss
un grafo e’ bipartito se e solo se non contiene cicli dispari
Accoppiamento di massima
cardinalita’ su grafo qualunque
G=(V,E), E = insieme delle coppie possibili
Max ijE xij:
ijS(i) xij  1  iN
xij{0,1}
b
e
d
a
S(i) e’ l’insieme degli archi incidenti nel nodo i
Ad esempio, la stella del nodo a e’
S(a) = {(a,b), (a,d), (a,c)}
g
f
c
Accoppiamenti e cicli dispari
• La soluzione intera ha al + k/2 archi  ciclo di k archi,
con k dispari
• La soluzione continua puo’ avere valore ½
su ogni arco del ciclo
Il vincolo
ijCk xij  k/2
per Ck ciclo dispari
e’ una disuguaglianza valida
Complessita’ del problema di matching
di massima cardinalita’ su grafi qualsiasi
Ck e’ un ciclo dispari di k archi
Il politopo
P={xij: ijS(i) xij  1 iN, ijCk xij  k/2 Ck, 0xij1}
e’ l’inviluppo convesso dei vettori di matching  risolvo il problema con la
PL se so descrivere tutti i vincoli..ma.. un numero esponenziale di cicli??
Se esiste una descrizione implicita polinomiale del politopo allora il problema
e’ in P anche se il politopo ha un numero esponenziale di facce:
Si puo’ dimostrare che esiste un algoritmo polinomiale per la generazione dei
tagli violati del tipo ijCk xij  k/2 (separazione)
 la separazione per i vincoli ijCk xij  k/2 e’ polinomiale
 il problema dell’accoppiamento di masssima cardinalita’ su grafi qualsiasi
e’ polinomiale
 Infatti esistono algortimi esatti ad hoc polinomiali
Un altro esempio:
il Minimum Spanning Tree
Dato
un grafo G=(V,E) non orientato
un vettore dei costi c:E→R+,
determinare il sottografo connesso
di costo minimo
Oss.
Un albero di copertura T e’ un sottografo
minimale per la connessione (rimuovendo un arco (i,j) qualsiasi da T l’insieme degli
archi rimanenti non forma + un grafo connesso, e si crea un TAGLIO, cioe’ una partizione
dei nodi V in due sottoinsiemi Q e V\Q con iQ, jV\Q non connessi fra loro),
massimale per l’assenza di cicli (aggiungendo a T un qualsiasi arco (i,j), il nuovo
grafo contiene un CICLO che include l’arco (i,j) )
Ciascuna delle due definizioni offre lo spunto per un modello valido di PLI.
I due modelli (per tagli e per cicli) sono entrambi corretti (definiscono lo stesso insieme di
punti interi) e con un numero esponenziale di vincoli, ma il loro rilassamento continuo (e il
bound associato) e’ molto diverso.
Un modello di PLI per il MST:
modello basato sui tagli
Min ijE cij xij tale che:
ijE xij = n-1
cardinalita’
iQ,jQ xij 1  QV, Q
connessione O(2n)
xij{0,1}
interezza
Data una soluzione intera x* (sottografo di n-1 archi) del rilassamento dato dalla cancellazione di
(alcuni dei) vincoli di connessione, esiste un algoritmo polinomiale per verificare la connessione
del sottografo:  visita sul grafo indotto dalla soluzione corrente x* .
Utilizzando i soli archi (i,j) con xij*=1, a partire dal nodo 1 visito tutti i nodi raggiungibili.
Se posso visitare tutti i nodi in V allora la soluzione e’ un albero di copertura.
Se al contrario x* viola il vincolo di connessione per un certo insieme di nodi Q, tale vincolo
viene introdotto esplicitamente nel modello, e si riottimizza il nuovo problema (procedimento
iterativo) fino alla soluzione del problema.
NB ad ogni iterazione risolvo un problema di PLI (potenzialmente richiede un B&B)…
Rilassamento continuo del
modello basato sui tagli per il MST
Min ijE cij xij tale che:
ijE xij = n-1
cardinalita’
iQ,jQ xij 1  QV
connessione O(2n)
xij[0,1]
rilassamento continuo dei vincoli di interezza
Supponiamo di risolvere il rilassamento continuo di un modello rilassato ottenuto esplicitando solo
un sottoinsieme dei vincoli di taglio.
E’ possibile generare eventuali vincoli violati per mezzo di una procedura di complessita’
polinomiale?
Data una soluzione del rilassamento continuo x* esiste una procedura polinomiale per verificare il
soddisfacimento di tutti i vincoli di connessione non presenti?
Premesso che
il problema di determinare la quantita’ massima di flusso che puo’ attraversare una rete
capacitata da un nodo s a un nodo t e’ un problema polinomiale (problema di flusso massimo)
• Utilizzo x* come capacita’ degli archi,
• Oriento la rete facendo corrispondere ad ogni [i,j]E una coppia di archi (i,j) e (j,i) A,
• Risolvo un problema di flusso massimo per ogni scelta della coppia di nodi s e t fra i nodi della
rete,
se ogni flusso massimo ha valore  1 allora TUTTI I vincoli di taglio sono rispettati, inclusi quelli non
esplicitati, altrimenti si genera il nuovo vincolo (e’ quello relativo al taglio di capacita’ minima
associato al flusso massimo <1) e si riottimizza
MST: modello basato sui tagli vs
modello basato sui sottocicli
Min ijE cij xij tale che:
ijE xij = n-1
cardinalita’
iQ,jQ xij 1  QV
connessione O(2n)
xij{0,1}
vincoli di interezza
Abbiamo appena visto che e’ possibile risolvere il rilassamento continuo di questa
formulazione in tempo polinomiale (poiche’ il problema di separazione lo e’).
Ho risolto il problema di MST? La soluzione ottenuta e’ intera?
Questo modello e’ una formulazione ideale dell’MST?
(l’ottimo del rilassamento continuo ha coordinate intere per ogni istanza?)
NO!!
Ma lo e’ una equivalente formulazione per la versione di PLI. Si tratta della formulazione basata
sui cicli in cui la connessione e’ imposta attraverso i vincoli di subtour breaking
Min ijE cij xij tale che:
ijE xij = n-1
cardinalita’
i,jQ xij  |Q|-1  QV
subtour breaking O(2n)
xij{0,1}
vincoli di interezza
MST: modello per tagli vs modello per sottocicli
Min ijE cij xij tale che:
ijE xij = n-1
iQ,jQ xij 1  QV
xij{0,1}
FORMULAZIONI
DI PLI
EQUIVALENTI
Min ijE cij xij tale che:
ijE xij = n-1
i,jQ xij  |Q|-1  QV
xij{0,1}
Pcut  Zn  Pcycle  Zn
Pcycle
Pcut
Min ijE cij xij tale che:
ijE xij = n-1
iQ,jQ xij 1  QV
xij[0,1]
Pcut
CONTIENE
Pcycle
Min ijE cij xij tale che:
ijE xij = n-1
i,jQ xij  |Q|-1  QV
xij[0,1]
Il politopo del rilassamento continuo della versione per tagli (detto Pcut)
CONTIENE (e’ una formulazione meno aderente)
il politopo del rilassamento continuo della versione per sottocicli (detto Pcycle).
Teorema:
Il politopo del rilassamento continuo della versione per tagli Pcut
CONTIENE il politopo del rilassamento continuo della versione per sottocicli Pcycle.
Min ijE cij xij t.c.:
ijE xij = n-1
iQ,jQ xij 1  QV
xij[0,1]
Pcut
CONTIENE Pcycle
Min ijE cij xij t.c.:
ijE xij = n-1
i,jQ xij  |Q|-1  Q  V
xij[0,1]
1) dimostriamo che una soluzione di Pcycle appartiene anche a Pcut
Consideriamo il vincolo di cardinalita’
ijE xij = n-1
Dato un qualunque sottoinsieme Q  V, gli archi in E si possono suddividere tra
quelli interni a Q, quelli interni a V\Q e quelli appartenenti al taglio d(Q)
Quindi ijE xij = i,jQ xij + i,jV/Q xij + i,jd(Q) xij = n-1
Sia x una soluzione in Pcycle,
allora x rispetta sia il vincolo di cardinalita’ che i vincoli di subtour breaking, per cui
i,jQ xij  |Q| - 1
i,jV/Q xij  |V\Q| -1
da cui i,jd(Q) xij  n-1 – (|Q|-1 + |V\Q|-1) = n -1 -n +2 =1
quindi x e’ anche una soluzione di Pcut
Teorema:
Il politopo del rilassamento continuo della versione per tagli Pcut
CONTIENE il politopo del rilassamento continuo della versione per sottocicli Pcycle.
2) dimostriamo che esiste una istanza per cui una soluzione di Pcut non appartiene a Pcycle
Consideriamo il grafo a sx in figura, in cui i costi degli archi sono rappresentati in blu.
L’albero ottimo ha costo 2 (es l’albero formato dagli archi {(ab), (bc), (bd), (ae)})
Nella figura a dx e’ rappresentata in rosso la soluzione ottima di Pcut, di costo 3/2.
Questo dimostra che Pcut non ha ottimo intero, quindi non e’ una formulazione ideale.
Inoltre la soluzione ottima di Pcut non e’ ammissibile per Pcycle, infatti la somma delle variabili
associate al ciclo C=(b, c, d) ha valore 5/2=2.5 > 2.
Non e’ possibile riallocare parte di tale valore sugli altri archi C senza aumentare il costo della
soluzione (gli archi del ciclo C sono gli unici a costo 0) ne’ ridurne il valore senza violare il vincolo
di cardinalita’ ( xijE = 4)
a
b
a
b
1
1
0
c
1/2
d
1
c
1/2
0
1
e
1/2
0
1
1/2
e
d
Teorema:
Pcycle e’ una formulazione ideale per MST:
Il politopo del rilassamento continuo della versione per sottocicli Pcycle
ha opttimo intero per ogni istanza di MST
Min ijE cij xij t.c.:
ijE xij = n-1
i,jQ xij  |Q|-1  Q  V
xij[0,1]
Pcycle
Domandiamoci se e’ facile risolvere Pcycle.
Una conseguenza importante della polimialita’ dell’algoritmo dell’elissoide, eguaglia
la complessita’ dell’algoritmo di separazione a quella del problema stesso
Quindi risolvere Pcycle e’ facile se, data x* soluzione ottima di un modello in cui solo un
sottoinsieme dei vincoli di subtour breaking e’ stato esplicitato, e’ facile determinare, se esiste,
un vincolo violato da x*.
Si tratta di cercare un insieme di nodi Q: jiQ x*ij |Q|
E’ sufficiente verificare che non esista un ciclo C la somma delle x* dei cui archi sia maggiore
di |C|-1.
Sia gij=1-x*ij 0 il costo dell’arco (ij). I vincoli di subtour elimination si possono riscrivere come
jiQ gij 1 e si risolva il problema di determinare per ogni nodo i il ciclo di costo minimo che passa
per i (problema polinomiale applicando Floyd-Warshall e leggendo I valori sulla diagonale
principale).
Se per ogni i tale valore e’ 1 allora x* soddisfa TUTTI I vincoli di subtour breaking.
esempio
• Il costo e’ proporzionale alla lunghezza
degli archi del disegno
• La soluzione iniziale
non e’ ammissibile, i
suoi archi non formano
2
un sottografo connesso
5
1
3
6
• Sul grafo rosso
effettuo una visita in
4
ampiezza a partire da
• I nodi marcati nella visita identificano
un nodo qualsiasi, es 1
quelli dell’insieme Q di un vincolo
violato relativo al taglio Q,V\Q
7
esempio
• A partire da 1 posso
raggiungere solo i nodi
1,2,3 e 4.
• Quindi esiste un taglio
formato dai due
insiemi di nodi,
N1={1,2,3,4} e
N2={5,6,7}
2
5
1
3
6
•
7
4
• Genero e aggiungo al modello il
vincolo x25 + x47  1 che e’ violato
dalla soluzione corrente
Risolvo nuovamente il
problema, e itero il
procedimento fino a
raggiungere l’
ammissibilita’, cioe’
un sottografo connesso