Programmazione Lineare
INTERA II
Scelta del bound
Rilassamento:
1. un problema la cui funzione obiettivo fornisce
una stima ottimistica di f(x) per ogni x
ammissibile di P
2. un problema la cui regione ammissibile include
F(P) (x eliminazione di vincoli)
3. 1 e 2
Un bound deve
• Essere veloce da calcolare (polinomiale)
• Fornire una stima del valore ottimo il piu’
possibile vicina (tight) al valore esatto
• Permettere di ottenere con poco sforzo
computazionale una soluzione euristica
ammissibile appartenente al sottoproblema, a
partire da quella non ammissibile del
rilassamento
Rilassamenti
•
Per eliminazione di vincoli. In particolare un rilassamento per eliminazione
di vincoli si dice rilassamento Combinatorio quando il problema rilassato e’
un problema di ottimizzazione combinatoria con una particolare struttura
che lo rende risolvibile in modo efficiente.
•
Convessificazione (si rende convesso un insieme discreto)
•
Surrogato (sostituisce a un set di vincoli la loro combinazione lineare)
•
Lagrangeano (elimina un set di vincoli, e aggiunge alla funzione obiettivo
una penalizzazione dipendente dalla violazione dei vincoli rimossi).
Introduce un nuovo problema relativo alla scelta delle penalita’ (quelle per
cui il rilassamento ottenuto fornisce il bound + tight) che si risolve
ottimizzando una funzione poliedrale (lineare a tratti in Rn). Questa tecnica
trova vasta applicazione nell’ottimizzazione non lineare.
Rilassamento per eliminazione di vincoli
un esempio di rilassamento combinatorio
Formulazione del TSP asimmetrico su
grafo orientato G=(N,A)
FS = Forward Star:
Sia xij {0,1}  arco (i,j)  A
Min ijA cijxij :
ijFS(i) xij =1  iN
jiBS(i) xji =1  iN
iS,jS xij 1  SN
xij {0,1}
{archi uscenti dal nodo i}
BS = Backward Star:
{archi entranti nel nodo i}
Vincoli di grado
Vincoli di connessione
Cosa resta?
Min ijA cijxij :
ijFS(i) xij =1  iN
jiBS(i) xji =1  iN
xij {0,1}  arco (i,j)  A
Problema di copertura per cicli di un grafo orientato
FACILE!!!
Infatti….. i vertici del politopo associato al rilassamento continuo
hanno coordinate intere. (Matrice totalmente unimodulare)
Questo e’ un esempio di un rilassamento COMBINATORIO
Rilassamento Surrogato
un esempio
Formulazione del TSP su grafo non-orientato G=(V,E)
Sia xij {0,1}  arco (i,j)  E, i<j
Min ijE cijxij :
ijS(i) xij =2  iV
iQ,jQ xij 1  QV
xij {0,1}
Vincoli di grado
Vincoli di connessione
S(i) = stella del nodo i = insieme degli archi incidenti nel nodo
Sostituisco ai vincoli di grado
la loro combinazione lineare a coefficienti 1
Min ijE cijxij :
Vincolo
di cardinalita’
ijS(i)
x
=
2

iV
n
2
ij
E
iQ,jQ xij 1  QV
xij {0,1}
Questo nuovo problema e’ facile (1-albero di costo minimo)
Rilassamento Lagrangeano
•
E’ un modo diverso per trarre vantaggio dalla struttura presente in alcuni vincoli
del modello (se fossero gli unici, il problema sarebbe “trattabile”)
•
Anziche’ eliminarli, vengono rimossi dal sistema dei vincoli e inseriti in funzione
obiettivo in modo tale da penalizzare la valutazione delle soluzioni non
ammissibili.
•
Si somma alla funzione obiettivo originaria il prodotto scalare tra un vettore di
penalita’ y (detti moltiplicatori Lagrangeani) e il vettore delle scarto associate ai
vincoli rimossi.
•
Esempio, sia P: max {cx : Axb, Exd, xZn} diventa, per una scelta del vettore
y0, RL(y): max {cx –y(b-Ax) : Exd, xZn}
•
Per ogni scelta di y0, RL(y) fornisce un upper bound di P
•
Si puo’ vedere il rilassamento per eliminazione di vincoli come un particolare
Lagrangeano che usa delle penalita’ nulle y=0.
Rilassamento Lagrangeano 2
• Si pone il problema della miglior scelta dei moltiplicatori y
• Per il problema di massimo, cerchiamo y tale RL(y) sia minimo (il +
possibile aderente al valore ottimo del problema di massimo originario).
• Si osservi che fissato x tale che (Exd, xZn), l’espressione cx –y(b-Ax) e’
lineare in y, quindi e’ un IPERPIANO NELLO SPAZIO DELLE y.
• Quindi RL(y) per y in Rm e’ una funzione descrivibile come il massimo al
variare di y di un fascio di iperpiani del tipo
{cx –y(b-Ax) : Exd, xZn}
Vedremo come ottimizzare
una funzione Poliedrale a tratti
quando parleremo di
NonLinearOoptimization (NLO)
Rilassamento continuo
(caso specifico di rilassamento per convessificazione)
L’insieme discreto delle soluzioni intere F(P) diventa un politopo
eliminazione dei vincoli di interezza
 x {0,1} diventa x  [0..1]
•
E’ per default il rilassamento implementato dai solver
•
Importante avere una buona formulazione del problema,
tale per cui non ci sia troppa differenza tra il politopo
che include F(P) (insieme delle soluzioni intere)
e l’inviluppo convesso di F(P) (Convex Hull, CH)
che corrisponde alla formulazione IDEALE di P
duality gap
•
Facilmente risolvibile (Programmazione Lineare)
•
Efficiente risoluzione di nodi figli/fratelli dell’albero attraverso la riottimizzazione
•
Tecniche sofisticate di generazioni di tagli (gestione centralizzata del pool dei tagli)
•
Permette di ottenere soluzioni approssimate a partire dalla soluzione frazionaria (non
sempre efficace con la programmazione 0/1)
Buona formulazione
E’ necessaria per avere un valore del
rilassamento continuo aderente a z*
Come si puo’ ottenere?
• Generazione di tagli
• Scelta dello spazio delle variabili decisionali: un
modello compatto e’ sempre preferibile? NO!!
P e’ una formulazione migliore di P
il bound fornito da P e’ piu’ aderente a z*
P e’ la formulazione IDEALE ( il bound e’ = z*)
Formulazione
IDEALE:
inviluppo convesso delle soluzioni intere:
Pre-processing:
Ci sono svariate tecniche per rafforzare
il modello in modo tale da ridurne le dimensioni,
eliminando vincoli ridondanti o viariabili
senza perdere le soluzioni ottime intere e migliorando il bound:
i solver + sofisticati le implementano automaticamente,
e permettono di scegliere se applicarle o meno,
e se applicarle solo al nodo radice o ai nodi di tutto l’albero
(trade off tra il rafforzamento del modello e il costo computazionale necessario)
I metodi + diffusi:
• Logical preprocessing
• Rafforzamento basato sui costi ridotti
• Introduzione di disuguaglianze valide
•
•
•
Variable Fixing basato sul valore dei costi ridotti delle
variabili nella soluzione ottima del rilassamento continuo
Sia x la soluzione ottima continua a un generico nodo dell’albero. Ci si domanda se sia
possibile fissare (per le soluzioni associate a tutto il sottoalbero di cui il nodo corrente e’ radice)
il valore di una variabile binaria che ha valore intero nel rilassamento continuo, in base al valore
del suo costo ridotto, in modo tale da non eliminare la soluzione ottima intera
La decisione e’ presa in base a un ragionamento”what if”.
Consideriamo un problema della forma max {cx: Ax=b, x0, xi {0,1} i}
e il suo rilassamento continuo
P: max {cx: Ax=b, xi 0, xi 1 i}
Considerando come vincoli espliciti anche i vincoli di segno, associando le variabili duali y ai
vincoli Ax=b, le variabili z ai vincoli x e (e e’ il vettore di tutti 1) e le variabili w ai vincoli x0 ,
il duale e’ il problema
D: min { yb+ze: yA+z–w=c, y,z,w 0}
•
Siano x*, e y*,z*,w* una coppia di soluzioni ottime per la coppia Primale Duale e definiamo
ci* = ci-y*Ai = zi*-wi* come il costo ridotto della variabile xi.
 il costo ridotto e’ la differenza dei prezzi ombra (var duali) dei vincoli di lower e upper bound
•
In base agli scarti complementari (primale di massimo)
–
–
–
•
•
se zi*>0 (quindi xi=1)  wi*=0
se wi*>0 (quindi xi=0)  zi*=0
se xi*<1 (quindi non e’ all’upper bund)  ci* 0 se xi*>0 (quindi non e’ al lower bound)  ci*0
se ci*>0  xi=1
se ci*<0  xi*=0
Sia z+ il valore dell’incumbent (un lower bound all’ottimo) e sia zLP il valore del rilassamento
continuo al nodo corrente. Sia xi*=0 (intera) e sia ci* = - wi*<0 (quindi wi*>0)
Se zLP + ci*  z+ allora xi puo’ essere FISSATA a 0.
Per esercizio, ripetere per il caso xi*=1
Logical preprocessing e probing
Savelsbergh Journal On Computing 6 (4) 1994, pp. 445-454
•
Si vuole ridurre la dimensione dei coefficienti nella matrice dei vincoli e dei
bounds sulle variabili, mantenendo inalterato l’insieme delle soluzioni intere.
Consideriamo il problema: Min {cx + hy: Ax+Gy≤b, x∊{0,1}n, y∊Rm , lj ≤ yj ≤ uj}
Indichiamo in ciascun vincolo con B+e B- (e rispettivamente C+ e C-) gli indici delle
variabili booleane (intere o continue) con coefficienti >0 e <0.
In tal modo, prendendo i valori assoluti dei coefficienti,
riformuliamo il vincolo iesimo aix + giy ≤ bi nella seguente forma:
Σ(j∊B+)ajxj - Σ(j∊B-)ajxj + Σ(j∊C+)gjyj - Σ(j∊C-)gjyj ≤ bi
•
Scopo: analizzare ciascun vincolo i per verificare se tale vincolo:
–
–
–
–
–
rende il problema inammissibile,
e’ ridondante,
puo’ restringere i bounds sulle variabili,
puo’ venir reso + aderente alla soluzione ottima intera modificandone i coefficienti,
forza il valore di una variabile nella soluzione ottima ad un unico valore 0 o 1.
Identificazione della non ammissibilita’
• Sia z il valore ottimo di Pi = min aix + giy :
Aix + Giy ≤ bi x∊{0,1}n, y∊Rm , lj ≤ yj ≤ uj
z e’ il minimo valore che assume la parte sx del vincolo i nella regione
ammissibile data dal sistema dei vincoli originari privato del vincolo i.
(indichiamo con Ai e Gi (bi) la matrice A e la matrice G (il vettore b) privati
dalla riga (componente) i-esima)
Se z > bi allora il vincolo i rende il problema inammissibile.
• Ovviamente il problema Pi e’ difficile quanto quello originario,
ma ogni lower bound su z che verifichi la condizione > bi , e’ sufficiente,
essendo zLB ≤ z, da cui zLB > bi  z > bi
Il lower bound + banale e + lasco e’ dato da:
zLB = - Σ(j∊B-)aj + Σ(j∊C+)gjlj - Σ(j∊C-)gjuj
in cui le variabili in B+ e C+ sono al lower bound,
e quelle in B- e C- sono all’upper bound
Identificazione della ridondanza
• Sia z il valore ottimo di Pi = max aix + giy :
Aix + Giy ≤ bi x∊{0,1}n, y∊Rm , lj ≤ yj ≤ uj
z e’ il massimo valore che assume la parte sx del vincolo i nella regione
ammissibile data dal sistema dei vincoli originario privato del vincolo i.
Se z ≤ bi allora il vincolo i e’ ridondante per il problema.
• Ovviamente il problema Pi e’ difficile quanto quello originario,
ma ogni upper bound su z e’ sufficiente,
essendo z ≤ zUB, da cui zUB ≤ bi  z ≤ bi
L’upper bound + banale e + lasco e’ dato da:
zUB = Σ(j∊B+)aj + Σ(j∊C+)gjuj - Σ(j∊C-)gjlj
in cui le variabili in B+ e C+ sono all’upper bound,
e quelle in B- e C- sono al lower bound
Miglioramento dei bounds
•
Sia zk il valore ottimo di Pki = min aix + Σ(jk, j∊C+) gjyj - Σ(j∊C-) gjyj
Aix + Giy ≤ bi x∊{0,1}n, y∊Rm , lj ≤ yj ≤ uj
zk e’ il minimo valore che assume la parte sx del vincolo i senza il contributo della
variabile yk, nella regione ammissibile data dal sistema dei vincoli originario
privato del vincolo i.
Ne segue che per soddisfare il vincolo i deve valere yk ≤ (bi-zk)/gik
Quindi l’upper bound su yk puo’ essere migliorato
se la nuova stima e’ inferiore alla precedente, cioe’ se (bi-zk)/gik≤uk
•
Ovviamente il problema Pki e’ difficile quanto quello originario,
ma ogni lower bound su zk e’ sufficiente,
essendo zkLB ≤ zk, da cui yk ≤ (bi- zkLB)/gik  yk ≤ (bi-zk)/gik
Il lower bound + banale e + lasco e’ dato da:
zkLB = - Σ(j∊B-)aj + Σ(jk,j∊C+)gjlj - Σ(j∊C-)gjuj
in cui le variabili in B- e C- sono all’upper bound,
e quelle in B+ e C+ sono al lower bound
esercizio
Fornire una procedura per il potenziale
miglioramento del lower bound lj su di una
variabile yj, secondo lo schema di
ragionamento descritto per il rinforzamento
dell’upper bound
Tecniche base di probing
(dedurre se il sistema di vincoli forza il valore di una variabile, caso xk=0)
•
•
Consideriamo un vincolo i e una variabile binaria xk in B+ e
sia zk il valore ottimo di Pki = min aix + giy :
Aix + Giy ≤ bi, xk=1, x∊{0,1}n, y∊Rm , lj ≤ yj ≤ uj
zk e’ il minimo valore che assume la parte sx del vincolo i nella regione ammissibile data dal
sistema dei vincoli originario privato del vincolo i
avendo forzato la variabile xk al valore 1.
Se zk > bi allora il problema Pk e’ inammissibile,
quindi xk=0 in ogni soluzione ammissibile al problema originario
e il suo valore puo’ essere fissato, riducendo cosi’ di 1 il numero delle variabili.
•
Ovviamente il problema Pk e’ difficile quanto quello originario,
ma ogni lower bound su zk e’ sufficiente,
essendo zLB ≤ zk, da cui zLB > bi  zk > bi
Il lower bound + banale e + lasco e’ dato da:
zLB = aik - Σ(j∊B-)aj + Σ(j∊C+)gjlj - Σ(j∊C-)gjuj
in cui le variabili in B+ e C+ sono al lower bound,
quelle in B- e C- sono all’upper bound, e xk=1.
Tecniche base di probing
(dedurre se il sistema di vincoli forza il valore di una variabile, caso xk=1)
•
•
Consideriamo un vincolo i e una variabile binaria xk in B- e
sia zk il valore ottimo di Pki = min aix + giy :
Aix + Giy ≤ bi, xk=0, x∊{0,1}n, y∊Rm , lj ≤ yj ≤ uj
zk e’ il minimo valore che assume la parte sx del vincolo i nella regione ammissibile data dal
sistema dei vincoli originario privato del vincolo i
avendo forzato la variabile xk al valore 0.
Se zk > bi allora il problema Pk e’ inammissibile,
quindi xk=1 in ogni soluzione ammissibile al problema originario
e il suo valore puo’ essere fissato, riducendo cosi’ di 1 il numero delle variabili.
•
Ovviamente il problema Pk e’ difficile quanto quello originario,
ma ogni lower bound su zk e’ sufficiente,
essendo zLB ≤ zk, da cui zLB > bi  zk > bi
Il lower bound + banale e + lasco e’ dato da:
zLB = - Σ(jk,j∊B-)aj + Σ(j∊C+)gjlj - Σ(j∊C-)gjuj
in cui le variabili in B+ e C+ sono al lower bound,
e quelle in B- e C- sono all’upper bound, e xk=0.
Un esempio numerico (L. Wolsey, Integer Programming cap. 7.6 pg103-105)
Consideriamo il seguente problema di PL (rilassamento continuo di un modello di PLI)
Il primo vincolo impone 5x1  15 +2x2 -8x3. In base ai bound si setta il valore delle variabili x2 e x3
nel modo + favorevole al vincolo (tale da max 15+2x2-8x3), da cui 5x1  15+2-8 = 9 da cui x1  9/5
che rende + stretto il bound precedente x1  3
max
2x1
+
x2
-
x3
tale
che
5x1
-
2x2
+
8x3

15
8x1
+
3x2
-
x3

9
x1
+
x2
+
x3

6
0

x1

3
0

x2

1
1

x3
Applicando lo stesso procedimento su x2 otteniamo
2x2  -15 + 5x1 + 8x3 e settiamo il valori delle variabili
in modo da minimizzare la parte dx del vincolo,
da cui 2x2  -15 + 5·0 + 8·1 = -7 ma questo non
migliora il bound attuale di x2 che e’ x2 0.
Lavorando su x3 abbiamo 8x3  15–5x1+2x2  15-0+2
da cui il nuovo bound su x3, x3  17/8
Ora si puo’ utilizzare il secondo vincolo (di ) per rendere + stringenti i lower bounds di x1 e x2
e l’upper bund di x3
Lavoriamo su x1: 8x1  9 -3x2 +x3  9 -3 +1 = 7 da cui x1  7/8 che rafforza x10
Invece, non e’ possibile migliorare i bounds di x2 e x3.
..poiche’ i bounds sono cambiati (sono diventati + stringenti), e’ opportuno reiterare il procedimento
Il primo vincolo mettendo in evidenza x3 e applicando i nuovi bounds su x2 e x1, produce
8x3  15 + 2·1 - 5·7/8 = 101/8 da cui x3101/64 che rafforza il precedente x317/8.
Continua…..
Ora il problema si presenta cosi’, gli stessi vincoli e funzione obiettivo ma i bounds sono + stretti.
Consideriamo il terzo vincolo, x1 +x2 +x3  6 e poniamo le variabili al loro upper bound.
Risulta 9/5 + 1 + 101/64 ~ 4.5 che e’ < 6, per cui ogni valore ammissibile delle variabili soddisfa il
terzo vincolo  il terzo vincolo e’ ridondante e si puo’ eliminare!!
max
2x1
+
x2
-
x3
tale
che
5x1
-
2x2
+
8x3

15
8x1
+
3x2
-
x3

9
x1
+
x2
+
x3

6
7/8

x1

9/5
0

x2

1
1

x3

101/64
Ora osserviamo i coefficienti di x2
Aumentandone il valore, ciascuno dei due vincoli
si allontana dal proprio right hand side
 al crescere di x2 i vincoli non vengono violati.
Poiche’ x2 ha coefficiente >0 in una funzione
obiettivo di massimo, nella soluzione ottima
x2 sara’ posto al suo upper bound x2=1
Quindi x2 puo’ essere fissata a tale valore
eliminandola dal modello.
x2=1!!
La variabile x3 ha coefficiente <0 nella funzione obiettivo di massimo,
Coefficiente >0 nella disuguaglianza di  e coefficiente <0
Il modello si riduce a
nella disuguaglianza di   in entrambi i casi diminuendo
Max {2x1 : 7/8  x1  9/5 }
il valore di x3 ci si avvicina alla (si resta nella) regione ammissibile
la cui soluzione ottima e’
 x3 va fissata al lower bound x3=1.
9/5, 1, 1.
x3=1!!
Disuguaglianze generate
dall’inferenza logica sulle variabili binarie
Consideriamo un sistema dato dalle seguenti 4 disuguaglianze, con x1..x4  {0,1}
Consideriamo il primo vincolo, 7x1 +3x2 -4x3 -2x4  1
Se x1=1 necessariamente sia x3 che x4 devono essere 1 
questo genera le implicazioni x1  x3 e x1  x4. Inoltre non e’ possibile che siano a 1 sia x1 che x2,
da cui segue il vincolo x1+x2  1.
Dal secondo vincolo segue che se x2=1 allora x1 deve essere 1, da cui l’implicazione x2  x1.
Inoltre x2 e x3 non possono essere contemporaneamente a 1, da cui x2+x3  1.
Dal terzo vincolo si deduce che o x4=1 oppure sia x2 che x3 sono 1, il che e’ imposto dal sistema:
x2 + x4  1, x3 + x4  1.
Dal quarto vincolo segue che x3  x1.
7x +
3x 4x
2x
1

1
-2x1 +
2
7x2
3
+
3x3
4
+
Vediamo cosa dedurre dai vincoli suddetti
-2x2 3x3
Il sistema x1+x2  1 e x2  x1 ha come unico
3x1
2x3
valore ammissibile per x2 il valore 0. x2=0!!
 possiamo fissare x2=0.
x4=1!!
x1 x2 x3 x4  {0,1}
Segue che x2+x4  1  x4=1
x3=x1!!
Da x1  x3 e x3  x1 segue x3=x1.
Sostituendo queste condizioni nel sistema, tutti i vincoli sono ridondanti,
Quindi sono ammissibili solo le soluzioni [1,0,1,1] e [0,0,0,1]
x4

6
6x4

-5

-1
Disuguaglianze valide
per problemi di PLI
Classi note di disuguaglianze
– Knapsack cover
– GUB cover
– Flow cover
– Cliques
– Implication
– Gomory
– Mixed-Integer rounding
Definizione di
Disuguaglianza Valida
Dato un problema di PLI del tipo
P = max cx: Axb, xZn,
il vincolo axa0 e’ una disuguaglianza valida per P
se e’ soddisfatto da ogni soluzione (intera) di P.
Dato x* punto di ottimo del rilassamento continuo di P:
PPL= max cx: Axb
una disuguaglianza valida si dice essere un taglio
se tale disuguaglianza è violata da x*
cioe’ si verifica che
ax*>a0
Tagli e bound
• Se aggiungo un taglio alla formulazione di un
problema di PLI, l’insieme delle soluzioni
ammissibili rimane tale, ma il bound associato al
rilassamento continuo delle nuova formulazione
migliora avendo eliminato il precedente punto di
ottimo del rilassamento continuo.
c
x*
Disuguaglianze valide generabili in
modo automatico su base algebrica
• Disuguaglianze di Chvatal
• Tagli di Gomory
sviluppo teorico anni 70  algoritmo
cutting plane
Inapplicate fino agli anni 90 (integrazione nei
metodi di B&B  Branch and Cut)
….dice Chvatal
dato il problema di programmazione lineare intera
(P) min{cx : Axb, xZ+n }
(A, b, c di componenti intere non negative)
ogni disuguaglianza del tipo ax  a0,
tale che a0 e ogni componente di a = [a1,..,an] sono
ottenute come parte intera inferiore della combinazione
conica delle righe di A|b con coefficienti uR+m
allora
la disuguaglianza ax  a0
è anche una disuguaglianza valida per P.
esempio
Interpretazione:
determinare un accoppiamento
di massima cardinalità
in un grafo completo di ordine 3 (K3)
a
x1
x2
c
b
x3
max x1 + x2 + x3
x1 + x2  1
x2 + x3  1
x1 + x3  1
x0, xiZ, i=1,…,3
Ogni variabile x rappresenta una coppia:
Il primo vincolo dice che il nodo a puo’ stare in coppia o con b o con c ma non con entrambi (o la
coppia (ab) o la coppia (ac) o nessuna delle due coppie).
Esistono 3 soluzioni ottime intere equivalenti corrispondenti a selezionare un qualsiasi arco
(coppia) (o x1 o x2 o x3 uguale a 1).
Quindi il valore dell’ottimo intero zI*=1 e xI*=[1,0,0] (ad esempio)
PERO’
L’ottimo del rilassamento continuo del modello e’ zLP*=3/2 con xLP*=[½, ½, ½].
L’ “integrality gap” (differenza relativa tra bound e ottimo) associato a questa istanza e’ (3/2-1)/1,
quindi dell’ordine del 50%
Vediamo come sia possibile rafforzare il modello in modo da migliorare (diminuire) il valore del
rilassamento continuo e restringere questo gap.
a
x1
x2
c
b
x3
max x1 + x2 + x3
x1 + x2
1
x2 + x3  1
x1
+ x3  1
x0, xiZ, i=1,…,3
u1
u2
u3
Prendiamo come moltiplicatori dei tre vincoli i tre numeri reali, u1=u2=u3=1/2,
e facciamo la combinazione conica dei vincoli, utilizzando come vettore dei coefficienti il vettore u
= [½, ½, ½]
NB ogni colore e’ associato a un vincolo in questo esempio
La disuguaglianza di Chvátal cosi’ ottenuta è:
½ + ½ x1 + ½ + ½ x2 + ½ + ½ x3  ½ + ½ + ½=1
cioè
x1 + x2 + x3  1.
Aggiungendo tale disuguaglianza al problema,
la soluzione ottima del rilassamento continuo e’ intera
e l’integrality gap e’ annullato
Un altro risultato teorico
Chvatal dimostra che:
•
e’ sufficiente un numero finito di tali disuguaglianze per definire la formulazione
ideale di P, attraverso vincoli che descrivono la convex hull (inviluppo convesso)
dei punti interi di P
(se avessimo la formulazione ideale del problema
si otterrebbe l’ottimo intero risolvendo un problema di PL!!)
•
Ma tali vincoli sono in numero esponenziale
•
dato x* ottimo del rilassamento continuo corrente,  sempre un vettore di
moltiplicatori u tale che la disuguaglianza uTAx≤uTb e’ violata da x*

QUINDI?
Poiche’ lo scopo e’ trovare l’ottimo intero z* e non ottenere l’intera descrizione
ideale di P, sono interessato alla generazione delle sole disuguaglianze che sono
anche tagli. Vale a dire, quelle disuguaglianze valide che sono violate dalla
soluzione ottima del rilassamento continuo corrente
Questa osservazione e’ alla base del metodo noto come Cutting Plane
Algoritmo Cutting Plane
Procedure metodo_piani_di_taglio:
begin
ottimo:= false;
repeat
x* := RisolviPL(P) ! x*= ottimo del rilassamento continuo
if x* intero then ottimo:= true
else begin
Genera_taglio (x*, a:a0);
Aggiungi_taglio(P, a:a0)
end
until ottimo
end.
Esiste una procedura per la generazione automatica
dei tagli di Chvatal, proposta da Gomory,
che sfrutta la soluzione di base associata a x*
Tagli di Gomory
Dato P= min {cx: Ax=b, x0 , xZn}, siano:
B la base ottima del rilassamento continuo,
A = [AB-1A] diventa la matrice [I, AN] nel tableau associato alla base ottima,
con AN=AB-1AN e sia b=AB-1b i valori della soluzione di base xB (xN=0).
Se x*=[xB,xN]=[b,0] non e’ intero  allora  un indice hB: x*h e’ frazionario ( bt e’ frazionario).
Sia t l’indice di riga del vincolo su cui xh ha il pivot (x*h= bt). Il vincolo nel tableau e’
xh + jN atj xj = bt
Passando ai troncamenti, genero un taglio per x* come segue, prendendo il floor dei coefficienti:
xh + jN atj xj ≤ bt 
Infatti, trattandosi di una disuguaglianza di tipo  ed essendo x0, vale sicuramente xh + jN atj xj
Il secondo step, che sostituisce a bt il suo floor, e’ dovuto ai vincoli di interezza.
Tale disuguaglianza e’ valida per le soluzioni intere di P, ma in x* e’ violata, perche’ vale
x*h + jN atj x*j > bt
per cui il nuovo vincolo e’ violato da x* essendo per hp x*h frazionario e x*j =0 per iN.
≤ bt .
Il nuovo vincolo viene aggiunto alla matrice corrente, in forma di uguaglianza, aggiungendo una apposita
variabile di slack intera, secondo una procedura detta di ORLATURA della matrice
che aggiunge una riga e una colonna ad ogni iterazione.
Queste operazioni, computazionalmente poco onerose nei solver di PLI ,sono implementate solo
nelle prime iterazioni poiche’ il procedimento perde via via di efficacia
(i tagli risultano progressivamente meno profondi).
Esempio numerico
min -x1 - x2
-3x1 + 12x2 + x3
= 30
6x1 - 3x2
+ x4 = 8
xiZ+
La soluzione ottima del problema continuo è xB = [62/21, 68/21], B = {1,2}, N={3,4} e
AB-1 = 1/21 1 4
2 1
Quindi il sistema originario, premoltiplicato per AB-1, diventa [I, AN] x = b
x1 + 1/21 x3 + 4/21 x4 = 62/21
x2 + 2/21 x3 + 1/21 x4 = 68/21
Applicando la procedura ad entrambe le componenti frazionarie della
soluzione, sia x1 che x2, produciamo rispettivamente per ciascuna riga,
i due nuovi vincoli corrispondenti ai tagli di Gomory :
x1 + 0 x3 + 0 x4  2,
x2 + 0 x3 + 0 x4  3
 x1 2 essendo 1/21= 4/21=2/21=0
 x2 3 essendo 62/21=2, 68/21=3
che sono rispettivamente violati dalla soluzione continua corrente x*=[62/21, 68/21, 0, 0]
esercizio
• Aggiungere i due tagli di gomory prodotti
all’esempio precedente, e risolvere
nuovamente il problema continuo con
Xpress
• Riapplicare la procedura fino
all’esaurimento (fino a quando non si
possono generare ulteriori tagli)
Aspetti pratici
•
Il solo algoritmo Cutting Plane utilizzando i tagli di Gomory non e’ efficace in
pratica: si osserva una progressiva perdita di efficacia dei tagli via via generati.
MA
diventa uno strumento molto potente se embedded in un Branch and Bound,
unito alla generazione di tagli ad hoc per particolari strutture del problema Branch
and Cut
IDEA alla base della procedura di Branch and Cut:
ad ogni nodo dell’albero di BB genero tagli fino a che sono efficaci. Se la
soluzione e’ ancora frazionaria allora Branch.
Vediamo alcuni tagli ad hoc per singoli problemi di ottimizzazione combinatoria, che
sono absati sulla specifica struttura del singolo problema.
Si osservi che il loro studio permette di applicare questi tagli a ogni formulazione di
PLI in cui un sottoinsieme dei vincoli abbia tale struttura.
Ad esempio molti modelli contengono almeno un vincolo del tipo axb.
Se la funzione obiettivo e’ di massimo e tutti i coefficienti sono non negativi, allora
posso applicare i tagli sviluppati per il problema di Knapsack intero
Tagli di tipo “cover”
riferimento per lo studio: folder “materiale didattico integrativo in pdf”
sul minisito del corso, file “cover-inequalities.pdf”
Dato un problema P a variabili binarie e un vincolo
del tipo axa0 di coefficienti ai0,
Si dice cover minimale del vincolo un gruppo di k
variabili Ck tali per cui:
- Se sono tutte a valore 1 il vincolo e’ violato
- Non appena una variabile assume valore 0, il
vincolo e’ soddisfatto
TH:
S(jCk) xj  k-1 e’ una disuguaglianza valida per P
Esempio
Max z = 9x1 + 5x2 + 6x3 + 4x4:
6x1 + 3x2 + 5x3 + 2x4  10
x3 + x4  1
-x1
+ x3
0
-x2
+ x4  0
xj {0,1}
Implementiamo il modello
con Xpress e risolviamone
il rilassamento continuo,
ottenendo x*
x*=[5/6,1,0,1] zL= 16,5
Sono cover derivate dal vincolo 1,
C3={1,2,4}, C2={1,3}
La cover C3 e’ violata da x* mentre C2 no. Quindi C2 NON SERVE!!!
Come generare solo cover violate? Cioe’ tali che
S(jCk) x*j > k-1
Risolvendo un nuovo problema di ottimizzazione!!!
Detto problema di SEPARAZIONE
Infatti il nuovo vincolo separa x* dalle soluzioni intere di P
model ModelName
uses "mmxprs";
OUTPUT
forward procedure stampa
declarations
x1, x2, x3, x4: mpvar
Objective:linctr
k : integer
end-declarations
6*x1 + 3*x2 + 5*x3 + 2*x4 <= 10
x3 + x4 <=1
-x1 + x3 <= 0
-x2 + x4 <= 0
x1 <= 1
x2 <= 1
x3 <= 1
x4 <= 1
Objective := 9*x1 + 5*x2 + 6*x3 + 4*x4
iterazione 1
16.5
x(1)=0.833333
x(2)=1
x(3)=0
x(4)=1
iterazione 2
16.5
x(1)=0.833333
x(2)=1
x(3)=0
x(4)=1
iterazione 3
14.5
x(1)=0.5 x(2)=1
x(3)=0.5 x(4)=0.5
k:=1
maximize(Objective)
stampa
x1 + x3 <= 1
maximize(Objective)
k:=k+1
stampa
x1 + x2 +x4 <=2
maximize(Objective)
k:=k+1
stampa
Aggiungere il vincolo relativo alla cover C2
non modifica la soluzione ottima continua x*
poiche’ C2 non e’ violata da x*
Al contrario, il vincolo relativo alla cover C1 e’ violato da x*
la nuova soluzione ottima del rilassamento continuo cambia,
e il bound associato migliora (da 16.5 scende a 14.5)
procedure stampa
writeln("\niterazione ",k,"\n", getobjval, "\tx(1)=", getsol(x1),
"\tx(2)=", getsol(x2), "\tx(3)=", getsol(x3), "\tx(4)=", getsol(x4))
end-procedure
end-model
Separazione per le cover inequalities:
Problema di ottimizzazione associato: come individuare tra le
tante cover, quella + violata da x* (ottimo continuo corrente).
Si chiama separazione perche’ il vincolo prodotto costituisce un
iperpiano di separazione tra x* e le soluzioni intere di P
Def.
C e’ una cover derivata dal vincolo ax  a0  SiC aixi* > a0
Usiamo una variabile binaria zi per denotare l’appartenenza
di un elemento i alla cover C.
Vogliamo
∑i aizi > a0
(una cover)
∑i (1-x*i) zi < 1
(il cui vincolo sia violato da x*)
NB: x* è nota quindi la uso come dato del nuovo problema,
inoltre vale (1-x*i)0 poiche’ 0  x 1.
Perche’???
Problema di Separazione per le cover:
..infatti, la cover vieta di avere contemporaneamente a 1 tutti i suoi
elementi, ed e’ espressa nella disuguaglianza
S (jCk) xj  k-1
Per cui e’ violata nel punto x* se S (jCk) xj* > k-1
Voglio imporre come vincolo, nel problema di ricerca della cover “migliore”,
che la cover sia violata da x*, cioe’ che per i valori correnti di x* valga
S(jCk) x*j > k-1
Opero delle trasformazioni di equivalenza a questa disuguaglianza:
- Sottraggo k a entrambi i membri:
- k + S(jCk) x*j > k - 1 - k = -1
- Porto –k dentro la sommatoria, distribuendo sui k addendi:
S(jCk) (x*j -1) > -1
- Moltiplico i membri per -1:
S(jCk) (1-x*j ) < 1
Giungendo cosi’ all’espressione finale del vincolo ∑i (1-x*i) zi < 1 che impone di
selezionare, attraverso la scelta delle variabili z che si sostituiscono
all’indicazione nella sommatoria (j Ck), solo cover violate da x*.
Poiche’ cerco la cover + violata, minimizzo l’espressione ∑i (1-x*i) zi
Cercare il vincolo di cover + violato (risolvere il problema di separazione),
equivale a risolvere il problema di ottimizzazione:
Min ∑i(1-x*i)zi tale che
∑i ai zi > a0, zi0,1 ristretto alle sole i con x*i>0,
Si tratta di un problema di zaino binario di dimensioni molto ridotte.
Nota che, poiche’ i coefficienti ai e a0 sono interi e le zi binarie,
il vincolo ∑i ai zi > a0 diventa ∑i ai zi  a0 +1
Ricordiamo che il problema di zaino binario,
anche se NP-Hard, e’ risolvibile in modo esatto da algoritmi pseudopolinomiali
(tramite la programmazione dinamica)
quindi e’ computazionalmente gestibile
risolvere la separazione con questa procedura.
Se il valore ottimo della funzione obiettivo è  1
allora TUTTI i vincoli di cover indotti dal vincolo [a,a0] sono soddisfatti da x*,
altrimenti
abbiamo generato la disuguaglianza di cover + violata,
la aggiungiamo all’insieme dei vincoli
e iteriamo la procedura risolvendo nuovamente il rilassamento continuo.
Lifting di un vincolo di cover (cenni)
•
Liftare una disuguaglianza significa renderla + efficacie, aumentando lo spazio in cui
tale vincolo e’ definito, in modo tale che sia sempre rispettato dalle soluzioni intere
del problema ma escluda ulteriori soluzioni frazionarie.
•
Un modo semplice di liftare una disuguaglianza originata da una cover C e’ quello di
aggiungere nella sommatoria degli elementi della parte sx della disuguaglianza, ogni
elemento il cui coefficiente “peso” sia maggiore di quello di tutti gli elementi della
cover stessa, lasciando inalterato il rhs.
•
Il principio su cui e’ basata questa operazione e’ il seguente: sia wmax il peso
dell’elemento piu’ pesante nella cover C. Allora, qualsiasi elemento di peso > wmax,
se sostituito a un qualunque elemento della cover C produce un insieme di elementi
il cui peso supera certamenete il limite a0 quindi il numero massimo di elementi
contemporaneamente selezionabili all’interno del nuovo insieme resta al piu’ k-1 (il
rhs). Non avremmo potuto affermare questa proprieta’ aggiungendo elementi di peso
inferiore a wmax.
•
•
•
Esempio: X={x∊{0,1}7 : 11x1+6x2+6x3+5x4+5x5+4x6+x7  19}
L’insieme di indici C={3,4,5,6} e’ una cover, poiche’ 6+5+5+4=20>19.
La disugaglianza x3+x4+x5+x6  3 puo’ essere rafforzata aggiungendo sia x1che x2. in
quanto il loro coefficiente e’ maggiore o uguale a quello degli elementi in C.
•
Il nuovo vincolo diventa x1+x2+x3+x4+x5+x6  3 quindi una disuguaglianza in R6
e si chiama “extended cover”
Lifting di un vincolo di cover (continua)
Una disuguaglianza di cover C, e’ un vincolo aderente al politopo (una faccia del
politopo) nello spazio RC (ristretto alle sole variabili della cover) se questa e’
minimale.
E’ possibile applicare il processo di lifting in modo da ottenere un vincolo in Rn tale che
resti il + possibile aderente al politopo?
Ordiniamo gli indici delle variabili  C secondo un criterio qualsiasi, e prendiamo la
prima (es. xh).
Il nostra scopo e’ liftare la cover ineq SIC xi  |C|-1 dallo spazio RC nello spazio RC{h} .
Dobbiamo rispondere alla seguente domanda: quale e’ il massimo valore che puo’
assumere il parametro b affinche’ il vincolo
b xh + SIC xi  |C|-1
sia una disuguaglianza valida per ogni soluzione intera del seguente problema di
zaino binario definito nello spazio RC{j} : ahxh + SIC aixi  Q?
Se xh=0 il vincolo e’ soddisfatto per ogni soluzione intera essendo C una cover.
Se xh=1, otteniamo il problema
Max b tale che b + Max {SIC xi tale che (ah + SIC ai xi)  Q }  |C|-1
Da cui b  |C|-1 - Max {SIC xi tale che SIC ai xi  Q-ah }
Il problema e’ un problema di zaino binario ma semplicissimo, risolvibile per ispezione,
data la particolare funzione obiettivo.
Lo si risolve, si setta il valore di b e si itera scegliendo la variabile successiva.
Lifting di un vincolo di cover (continua)
Esempio: sia P il problema di zaino binario
P={max wx tale che x∊{0,1}7 e 11x1+6x2+6x3+5x4+5x5+4x6+x7  19}
e sia C={x3, x4, x5, x6} una cover per P  SiCxi  |C|-1  x3+x4+x5+x6  3 e’
una disuguaglianza valida per P ed e’ “facet defining” per il problema di
zaino in RC dato da PC={max wx, tale che x∊{0,1}4 e 6x3+5x4+5x5+4x6 19}
Cio’ significa che il vincolo e’ aderente alla frontiera di PCR4 per un
sottospazio di dimensione 4-1=3, e quindi ne definisce una “faccetta”.
Partiamo dalla variabile x1. Vogliamo liftare il vincolo in RC{1} secondo la
procedura illustrata in precedenza.
Occorre determinare b*= max b: b + x3+x4+x5+x6  3 e’ soddisfatta da ogni
soluzione intera del problema 6x3+5x4+5x5+4x6 19-11=8, con x∊{0,1}4
Il valore della soluzione ottima di max {x3+x4+x5+x6 tale che 6x3+5x4+5x5+4x6
8, con x∊{0,1}4 } e’ 1 (si verifica facilmente per ispezione), quindi b*= max
b: (b + 1)  3 vale 2, e otteniamo il vincolo di cover liftato
2x1 +x3+x4+x5+x6  3
Lifting di un vincolo di cover (continua)
Ora liftiamo il vincolo 2x1 +x3+x4+x5+x6  3 rispetto alla variabile x2.
Ripetendo I passaggi precedenti, dobbiamo risolvere il problema di zaino
max {2x1+x3+x4+x5+x6 tale che 11x1+6x3+5x4+5x5+4x6  19-6 = 13, xi{0,1} }
Ci sono diverse soluzioni intere con valore ottimo 2, quali:
x1=1, x3=x4=x5=x6=0,
x1=x3=x6=0, x4=x5=1,
x1=x4=x5=0, x3=x6=1,
Quindi b*= max b: b + 2 3  b*= 1, il che produce la seguente cover inequality liftata in RC{1,2}
2x1+x2+x3+x4+x5+x6  3.
Iteriamo il procedimento su x7 per liftare il vincolo in R7.
Il problema max {2x1+x2+x3+x4+x5+x6 t.c. 11x1+6x2+6x3+5x4+5x5+4x6  19-1=18, xi{0,1} }
ha come valore ottimo 3, associato a diverse soluzioni intere quali
x2=x3=x4=1, x1=x5=x6=0, ma anche x4=x5=x6=1, x1=x2=x3=0.
Quindi b*= max b: b + 3 3  b*= 0, per cui il coefficiente massimo che puo’ assumere x7 nella
cover liftata e’ 0 per cui il vincolo di cover mantiene la stessa formulazione anche in R 7.
Il procedimento ha termine con la garanzia che il vincolo ottenuto e’ aderente a CH(X(P)), la
combinazione convessa di tutte le soluzioni intere del problema di zaino binario di partenza.
Ogni diverso ordinamento delle variabili C da potenzialmente origine a un vincolo differente.
Quindi il procedimento andrebbe iterato per ogni possibile ordinamento. Svolgere questa parte
per esercizio.
Riferimento bibliografico, L. wolsey, Integer Programming, cap 9.3, pg147-151
esercizio
• Utilizzando Xpress, aggiungere la disuguaglianza di cover violata
all’esempio numerico HL9.6 e rivalutare il valore del rilassamento
continuo.
• Ci sono altre cover violate? Fornire una dimostrazione della risposta
• Liftare ciascuna delle cover ottenute, secondo tutti gli ordinamenti
delle variabili non appartenenti alla cover
• Implementare in Xpress un procedimento iterativo che dato un tale
ordinamento, implementa il lifting.
Un altro esempio di taglio generabile automaticamente dai solver:
i tagli di tipo clique
Esercizio: Formula il problema
di Air Crew Scheduling settimanale per il seguente input:
Voli da a
1
2
3
4
5
6
7
A
B
C
C
D
A
B
B
C
D
A
A
B
C
td
ta
8:00
10:00
13:00
15:00
15:00
17:00
11:00
9:00
11:00
14:00
16:00
16:00
18:00
12:00
freq (settimanale)
1,2,3,4,5
1,2,3,4,5,6,7
1,2,3,4,5,6
1,2,3,4,5,6,7
1,2,3,4,5,7
1,2,3,4,5,6,7
1,2,3,4,5
Regole:
•
1h fra 2 voli consecutivi;
•
non + di 6h di volo al giorno in un turno;
•
A,C,D sono basi,
•
turni di al + 2 giorni.
Sono (alcuni fra i) possibili turni (colonne del modello di Set Partitioning) le
seguenti sequenze di voli, indicati come
(T_numturno, giorno inizio) = giorno: volo1,volo2.. / giorno+1: volo1,volo2..
•
•
•
•
•
•
•
•
(T1,i) = i: v1,v7,v4
(T2,i) = i: v6 / i+1: v2, v3, v5
(T3,i) = i: v1, v2, v3, v5
(T4,i) = i: v6 / i+1: v7,v4
(T5,i) = i: v1, v2, v4
(T6,i) = i: v6 / i+1: v7, v3, v5
(T7,i) = i: v5, v6 / i+1: v2, v3
(T8,i) = i: v6 / i+1: v2, v4
per i=1,2,3,4,5
per i=1,2,3,4
per i=1,2,3,4,5
per i=7,1,2,3,4
per i=1,2,3,4,5
per i=7,1,2,3,4
per i=7
per i=1,2,3,4,5,6,7
OSSERVAZIONE:
Ci sono delle incompatibilita’ tra coppie di turni in giorni specifici,
Esempio: non possono essere nella stessa soluzione (T1,2) e (T4,1) perche’
entrambi coprono i voli v7 e v4 del giorno 2.
Conclusioni
•
Se disponiamo di una formulazione ideale di un problema di PLI
I
P = min {cx: xX} = min {cx: Axb, x0, x intero} = min {cx : x PPLZn}
(cioe’ conosciamo l’espressione algebrica esplicita dei vincoli del politopo corrispondente all’inviluppo convesso
delle soluzioni intere P*PPL ),
Potremmo risolvere PI con l’algoritmo del simplesso applicato a P*.
Quando tale formulazione presenta un numero esponenziale di vincoli tale procedimento e’ polinomiale SE E
SOLO SE
E’ possibile risolvere in tempo polinomiale il problema di separazione:
dato x* ottimo di min {cx: Axb, x0} = min {min {cx : x PPL} allora
o x*Zn e il problema PI e’ risolto all’ottimo
oppure determino un taglio ax  a0 tale che ax  a0  xX ma ax*>a0
Ad esempio
La formulazione per subtour breaking del problems di MST e’ una formulazione ideale, e i vincoli di subtour
breaking sono separabili con complessita’ polinomiale,
Ne consegue che il problema di MST appartiene alla classe P (problemi “facili”)
In caso contrario, se solo alcune delle famiglie esponenziali di vincoli che defiscono PI sono facilmente separabili,
nell’ambito di un processo di Branch & Bound (B&B) e’ opportuno generare tutti quei vincoli che tagliano la
soluzione continua ottima attraverso la soluzione del problema di separazione, e solo in seguito ricorrere al
branching, avendo cosi’ un bound di migliore qualita’ e riducendo i tempi del B&B
Tale processo ibridizza il Branch and Bound con la generazione di tagli tipica dell’algoritmo di Cutting Plane, e
prende il nome di Branch & Cut (B&C)
Nel Branch and Cut alcuni tagli possono essere validi solo ocalmente per il nodo corrente, poiche’ sono
condizionati dagli specifici branching associati al cammino dal nodo rdice al nodo corrente.
Altri tagli, di carattere globale sono validi in ogni ndoo dell’albero.
La gestione efficiente di questi tagli richiede algoritimi ad hoc per non appesantire il solver.