soluzione

laboratorio 2
Ottimizzazione
Prof. E. Amaldi
1) Risoluzione del rilassamento continuo del problema del commesso viaggiatore
Sia G = (V, E) un grafo orientato completo, con un costo cij ∈ R associato a ciascun arco (i, j) ∈ E. Si consideri la seguente formulazione (Dantzig, Fulkerson, Johnson, 1959) del
problema del commesso viaggiatore asimmetrico (Asymmetric Traveling Salesman Problem,
ATSP)
X
cij xij
(1)
(DF J)
min
(i,j)∈E
X
s.t.
xij = 1
i∈V
(2)
xji = 1
i∈V
(3)
xij ≥ 1
S ⊂ V, 1 ∈ S, |S| > 1
(4)
(i, j) ∈ E,
(5)
(i,j)∈δ + (i)
X
(i,j)∈δ − (i)
X
(i,j)∈δ + (S)
xij ∈ {0, 1}
dove i vincoli (4) rappresentano le cut-set inequalities CUT, e il suo rilassamento continuo
X
(DF J 0 )
min
cij xij
(6)
(i,j)∈E
s.t.
X
xij = 1
i∈V
(7)
xji = 1
i∈V
(8)
xij ≥ 1
S ⊂ V, 1 ∈ S, |S| > 1
(9)
(i,j)∈δ + (i)
X
(i,j)∈δ − (i)
X
(i,j)∈δ + (S)
xij ≥ 0
(i, j) ∈ E.
(10)
a) Si esprimano le formulazioni in linguaggio AMPL.
b) Si osservi sperimentalmente che, per l’istanza proposta, il rilassamento continuo DF J 0
fornisce soluzioni intere.
c) Dato che i vincoli (9) sono in numero esponenziale nel numero di nodi |V |, si implementi
in AMPL un metodo dei piani di taglio per risolvere DF J 0 , in cui i vincoli (9) vengono
introdotti iterativamente solo se violati.
Suggerimento: data una soluzione frazionaria di (7)-(10), con solamente un sottoinsieme
di disuguaglianze (9) generate finora, occorre individuare, se esiste, una diseguaglianza
violata appartenente alla classe (9) (se non esistono vincoli violati, la soluzione `e ammissibile e quindi ottimale). Si riconduca il problema di separazione alla risoluzione di una
sequenza di |V | − 1 problemi di taglio di capacit`a minima in un grafo con capacit`a e nodi
sorgente e destinazione scelti in modo opportuno.
Documento preparato da S. Coniglio, L. Taccari, V. Danese
1
laboratorio 2
Prof. E. Amaldi
Ottimizzazione
2) Formulazione alternativa
Si consideri la seguente formulazione estesa (Miller, Tucker, Zemlin) per il problema del
commesso viaggiatore asimmetrico:
X
cij xij
(11)
(MTZ) min
(i,j)∈E
s.t.
X
xij = 1
∀j ∈ V
(12)
xij = 1
∀i ∈ V
(13)
∀ (i, j) ∈ E, i 6= 1, j 6= 1
(14)
∀ (i, j) ∈ E
(15)
i∈V
X
j∈V
ui − uj ≤ (|V | − 1) (1 − xij ) − 1
xij ∈ {0, 1}
Quanti vincoli contiene? Come `e rispetto a (DF J)?
d) Si interpreti il significato delle variabili u.
e) Si esprima la formulazione (M T Z) in AMPL e si osservi la soluzione del rilassamento
continuo (M T Z 0 ). E’ ancora intera?
Traccia della soluzione
Segue il file k6 dir.dat per un grafo completo con 6 nodi (valore ottimo della funzione
obiettivo: 135):
#k6_dir.dat
data;
set V := a b c d e f;
param c:
a
b
c
d
e
f
;
a
.
02
04
99
87
54
b
c
d e
02 04 99
.
03 68
03
. 99
68 99
.
74 98 02
54 84 03
f
87
74
98
02
.
04
:=
54
54
84
03
04
.
Segue il file dfj-i.mod contenente il modello parziale della formulazione (DFJ), da com-
Documento preparato da S. Coniglio, L. Taccari, V. Danese
2
laboratorio 2
Ottimizzazione
Prof. E. Amaldi
pletare:
#dfj-i.mod
set V ordered;
param n := card {V};
set POWERSET := 0 .. (2**n - 1);
set S {k in POWERSET} := {i in V: (k div 2**(ord(i)-1)) mod 2 = 1};
set E := {i in V, j in V: ord(i) <> ord(j)};
#add parameters declarations
#add variables declarations
#add model
Osserviamo che un generico sottoinsieme S ⊆ V pu`o essere rappresentato con una stringa di
n = |V | caratteri a valori in {0, 1}, in cui al posto i compare 1 se il nodo i `e appartenente ad S, 0
altrimenti. La rappresentazione in base 2 dei numeri interi compresi tra 0 e 2n −1 fornisce quindi
tutte le possibili stringhe di n caratteri a valori in {0, 1}. Pertanto, per costruire tutti i possibili
sottoinsiemi di V, col comando set POWERSET := 0 .. (2**n - 1); definiamo l’insieme degli
indici da 0 a 2n − 1. Scegliamo di rappresentare il sottoinsieme S [k] con la codifica del numero
k in base 2. S [k] sar`
a quindi formato dai nodi la cui cifra corrispondente `e pari a 1. Poich`e in
base 2 la cifra i-esima `e pari a 1 se e solo se
(k div 2**(ord(i)-1)) mod 2 = 1
allora
set S {k in POWERSET} := {i in V: (k div 2**(ord(i)-1)) mod 2 = 1};
Il modello viene risolto lanciando il file dfj.run:
model dfj.mod
data k6_dir.dat
option solver cplexamp;
solve;
display x;
printf "Rilassamento continuo\n\n";
option relax_integrality 1;
solve;
display x;
Documento preparato da S. Coniglio, L. Taccari, V. Danese
3
laboratorio 2
Ottimizzazione
Prof. E. Amaldi
Segue un’implementazione parziale del metodo dei piani di taglio, nel file cuttingplanes-i.run:
#cuttingplanes-i.run
model primal.mod;
model separation.mod;
data k6_dir.dat;
option solver cplex;
problem primal:
#list vars, objective function, constraints names
problem separation:
#list vars, objective function, constraints names
let nc := 0;
repeat {
printf "Iteration %s, solving current relaxation\n", nc;
solve primal;
display x;
let s := first(V);
for {tt in V: tt <> s} {
let t := tt;
printf "Generating cut for s=%s, t=%s\n", s, t;
solve separation;
if (separationObj < 1) then {
let nc := nc + 1;
#add the cut to the constraints of the primal problem
printf "Added the cut: "; display CUT[nc];
break;
}
}
} while (separationObj < 1);
solve primal;
display x;
printf "Cuts added = %s\n", nc;
Segue il file MTZ-i.mod contenente il modello parziale della formulazione (MTZ), da com-
Documento preparato da S. Coniglio, L. Taccari, V. Danese
4
laboratorio 2
Ottimizzazione
Prof. E. Amaldi
pletare:
# INSIEMI
set V ordered;
set V_ridotto := {i in V : ord(i) > 1};
param n := card(V);
set E := {i in V, j in V: ord(i) <> ord(j)};
# PARAMETRI
# VARIABILI
# FUNZIONE OBIETTIVO
# VINCOLI
Il modello viene risolto lanciando il file MTZ.run:
model MTZ.mod
data k6_dir.dat
option solver cplexamp;
solve;
display x;
display u;
option relax_integrality 1;
solve;
display x;
display u;
Documento preparato da S. Coniglio, L. Taccari, V. Danese
5
laboratorio 2
Ottimizzazione
Prof. E. Amaldi
Soluzione
a) Segue il file di modello dfj.mod completo:
# INSIEMI
set V ordered;
param n := card(V);
set POWERSET := 0 .. (2**n -1);
set S{k in POWERSET} := {i in V : (k div 2**(ord(i)-1)) mod 2 = 1};
set E := {i in V, j in V: ord(i) <> ord(j)};
# PARAMETRI
param c{E}, >=0;
# VARIABILI
var x{E}, binary;
# FUNZIONE OBIETTIVO
minimize costi:
sum{(i,j) in E} c[i,j]*x[i,j];
# VINCOLI
subject to grado_in{i in V}:
sum{j in V: (j,i) in E} x[j,i] = 1;
subject to grado_out{i in V}:
sum{j in V: (i,j) in E} x[i,j] = 1;
subject to tagli{k in POWERSET diff {2**n-1}: (k div 2**(1-1)) mod 2 = 1
and card(S[k]) > 1}:
sum{i in S[k], j in V diff S[k]: (i,j) in E} x[i,j] >= 1;
La soluzione `e la seguente:
Documento preparato da S. Coniglio, L. Taccari, V. Danese
6
laboratorio 2
Prof. E. Amaldi
Ottimizzazione
CPLEX 11.2.0: optimal integer solution; objective 135
17 MIP simplex iterations
0 branch-and-bound nodes
x [*,*]
:
a
b
c
d
e
f
:=
a
.
0
0
0
0
1
b
0
.
1
0
0
0
c
1
0
.
0
0
0
d
0
1
0
.
0
0
e
0
0
0
1
.
0
f
0
0
0
0
1
.
;
Rilassamento continuo
CPLEX 11.2.0: optimal solution; objective 135
30 dual simplex iterations (0 in phase I)
x [*,*]
:
a
b
c
d
e
f
:=
a
.
0
0
0
0
1
b
0
.
1
0
0
0
c
1
0
.
0
0
0
d
0
1
0
.
0
0
e
0
0
0
1
.
0
f
0
0
0
0
1
.
;
Come si pu`
o notare, per l’istanza considerata il rilassamento continuo (DF J 0 ) fornisce
ancora una soluzione intera.
c) L’algoritmo di generazione di piani di taglio procede come segue:
1) Si risolva DF J 0 , dove tutte le cut-set inequalities sono rilassate, non venendo inserite
nella formulazione. Sia x∗ una sua soluzione ottima.
2) Si trovi una disuguaglianza violata, cercando un taglio nel grafo i cui lati abbiano
` dunque sufficiente cercare un taglio a
peso complessivo strettamente inferiore a 1. E
peso minimo e verificare se il suo valore `e inferiore a 1. Si osservi che questo problema
`e risolvibile in tempo polinomiale se i pesi del grafo sono nonnegativi. Troviamo un
taglio di peso minimo cercando un s-t taglio (un taglio che separa s da t) n volte, una
per ogni possibile vertice t ∈ V fissato un s qualsiasi, risolvendo il seguente problema
di programmazione lineare:
X
(CU T − SET )
min
x∗ij αij
(i,j)∈E
s.t. πs − πt ≥ 1
αij ≥ πi − πj
(i, j) ∈ E
αij ≥ 0
(i, j) ∈ E,
dove α `e il vettore di incidenza di un taglio del grafo, e π `e il vettore dei cosiddetti
potenziali di nodo. Il potenziale di s deve essere almeno di un’unit`a maggiore rispetto
al potenziale di t, e quindi, in un grafo completo, il vettore α non pu`o contenere solo
` facile mostrare che in ogni soluzione ottima del problema, α `e un vettore binario,
0. E
ed in particolare i componenti del vettore α a valore 1 indicano i lati appartenenti
al taglio. Se il valore del taglio minimo `e inferiore a 1, i lati appartenenti al taglio
violano la cut-set inequality.
Documento preparato da S. Coniglio, L. Taccari, V. Danese
7
laboratorio 2
Ottimizzazione
Prof. E. Amaldi
Si noti che la formulazione corrisponde al duale del problema di flusso massimo tra
un nodo s e un nodo t.
(M AX − F LOW )
max φ

if i=s

φ
X
X
fji = −φ if i=t i ∈ V
fij −
s.t.


(j,i)∈δ − (i)
(i,j)∈δ + (i)
0
else
∗
fij ≤ xij
(i, j) ∈ E
3) Si aggiunga la disuguaglianza alla formulazione e la si risolva nuovamente.
4) Si iteri fin che il valore del taglio di capacit`a massima non risulti non inferiore a 1.
File di modello primal.mod, predisposto per l’inserimento di nuove disuguaglianze:
set V ordered;
param n := card{V};
set E := {i in V, j in V: ord(i) <> ord(j)};
param nc >=0;
set CUT{1..nc} within V;
param c{E}, >=0;
var x{E}, >=0, <=1;
minimize obj:
sum{(i,j) in E} c[i,j]*x[i,j];
subject to grado_in{i in V}:
sum{j in V: (j,i) in E} x[j,i] = 1;
subject to grado_out{i in V}:
sum{j in V: (i,j) in E} x[i,j] = 1;
subject to cutSetInequalities{e in 1..nc}:
sum{i in CUT[e], j in V diff CUT[e]: (i,j) in E} x[i,j] >= 1;
Documento preparato da S. Coniglio, L. Taccari, V. Danese
8
laboratorio 2
Ottimizzazione
Prof. E. Amaldi
Il problema di separazione `e definito nel file separation.mod:
param s symbolic in V, default 1;
param t symbolic in V;
var alpha{E}, >=0;
var pi{V};
minimize separationObj:
sum{(i,j) in E} (x[i,j]*alpha[i,j]);
subject to rootPotential:
pi[s] - pi[t] >= 1;
subject to potentials{(i,j) in E}:
alpha[i,j] >= pi[i] - pi[j];
Segue il file cuttingplanes.run, con l’implementazione completa del metodo dei piani di
taglio:
#cuttingplanes.run
model primal.mod;
model separation.mod;
data k6_dir.dat;
option solver cplexamp;
problem primal:
x, obj, grado_in, grado_out, cutSetInequalities;
problem separation:
alpha, pi, separationObj, rootPotential, potentials;
let nc := 0;
repeat {
printf "Iteration %s, solving current relaxation\n", nc;
solve primal;
display x;
let s := first(V);
for {tt in V: tt <> s} {
let t := tt;
printf "Generating cut for s=%s, t=%s\n", s, t;
solve separation;
if (separationObj < 1) then {
let nc := nc + 1;
let CUT[nc] := {i in V: pi[i] = pi[s]};
printf "Added the cut: "; display CUT[nc];
break;
}
}
} while (separationObj < 1);
solve primal;
display x;
printf "Cuts added = %s\n", nc;
Documento preparato da S. Coniglio, L. Taccari, V. Danese
9
laboratorio 2
Ottimizzazione
Prof. E. Amaldi
La soluzione ottima trovata col metodo `e:
Documento preparato da S. Coniglio, L. Taccari, V. Danese
10
laboratorio 2
Ottimizzazione
Prof. E. Amaldi
Iteration 0, solving current relaxation
CPLEX 11.2.0: optimal solution; objective 18
10 dual simplex iterations (0 in phase I)
x [*,*]
:
a
b
c
d
e
f
:=
a
.
1
0
0
0
0
b
0
.
1
0
0
0
c
1
0
.
0
0
0
d
0
0
0
.
1
0
e
0
0
0
0
.
1
f
0
0
0
1
0
.
;
Added the cut: set
CUT[nc] := a b c;
Iteration 1, solving current relaxation
CPLEX 11.2.0: optimal solution; objective 118
2 dual simplex iterations (0 in phase I)
x [*,*]
:
a
b
c
d
e
f
:=
a
.
0
0
0
0
1
b
0
.
1
0
0
0
c
0
1
.
0
0
0
d
0
0
0
.
1
0
e
0
0
0
1
.
0
f
1
0
0
0
0
.
;
Added the cut: set
CUT[nc] := a d e f;
Iteration 2, solving current relaxation
CPLEX 11.2.0: optimal solution; objective 119
6 dual simplex iterations (4 in phase I)
x [*,*]
:
a
b
c
d
e
f
:=
a
.
0
1
0
0
0
b
0
.
0
0
0
1
c
0
1
.
0
0
0
d
0
0
0
.
1
0
e
0
0
0
1
.
0
f
1
0
0
0
0
.
;
Added the cut: set
CUT[nc] := a b c f;
Iteration 3, solving current relaxation
CPLEX 11.2.0: optimal solution; objective 135
9 dual simplex iterations (4 in phase I)
x [*,*]
:
a
b
c
d
e
f
:=
a
.
0
1
0
0
0
b
0
.
0
1
0
0
c
0
1
.
0
0
0
d
0
0
0
.
1
0
e
0
0
0
0
.
1
f
1
0
0
0
0
.
;
Cuts added = 3
Documento preparato da S. Coniglio, L. Taccari, V. Danese
11
laboratorio 2
Ottimizzazione
Prof. E. Amaldi
La soluzione `e stata trovata generando solo 3 disuguaglianze, in luogo di 2|V |−1 − 1 = 31.
d) In questa formulazione estesa viene aggiunto un numero lineare di variabili (n = |V |) e
un numero polinomiale di vincoli (nel caso di un grafo orientato completo (n − 1) (n − 2)
vincoli, in generale O n2 ). Pertanto `e una formulazione compatta.
Le variabili {ui }i∈V possono essere interpretate come l’ordine di visita dei nodi del grafo.
Segue il file di modello MTZ.mod completo:
# INSIEMI
set V ordered;
set V_ridotto := {i in V : ord(i) > 1};
param n := card(V);
set E := {i in V, j in V: ord(i) <> ord(j)};
# PARAMETRI
param c{E}, >=0;
# VARIABILI
var x{E}, binary;
var u{V};
# FUNZIONE OBIETTIVO
minimize costi:
sum{(i,j) in E} c[i,j]*x[i,j];
# VINCOLI
subject to grado_in{i in V}:
sum{j in V: (i,j) in E} x[i,j] = 1;
subject to grado_out{i in V}:
sum{j in V: (j,i) in E} x[j,i] = 1;
subject to ordine_visita{i in V_ridotto, j in V_ridotto:(i,j) in E}:
u[i] - u[j] <= (n-1)*(1-x[i,j]) -1;
###########
# NB: questi due vincoli non sono necessari per la correttezza della formulazione
# (sono solo per avere la corrispondenza tra la variabile u[i] e
# l’ordine di visita del nodo i)
subject to nodo_1:
u[first(V)] = 1;
subject to altri_nodi{i in V: i <> first(V)}:
u[i] >= 2;
############
Confrontiamo la soluzione fornita dal rilassamento continuo con il problema di partenza:
Documento preparato da S. Coniglio, L. Taccari, V. Danese
12
laboratorio 2
Ottimizzazione
Prof. E. Amaldi
CPLEX 11.2.0: optimal integer solution; objective 135
50 MIP simplex iterations
10 branch-and-bound nodes
8 implied-bound cuts
1 mixed-integer rounding cut
x [*,*]
:
a
b
c
d
e
f
:=
a
.
0
1
0
0
0
b
0
.
0
1
0
0
c
0
1
.
0
0
0
d
0
0
0
.
1
0
e
0
0
0
0
.
1
f
1
0
0
0
0
.
;
u [*] :=
a 1
b 3
c 2
d 4
e 5
f 6
;
CPLEX 11.2.0: optimal solution; objective 18
22 dual simplex iterations (0 in phase I)
x [*,*]
:
a
b
c
d
e
f
:=
a
.
0.8
0.2
0
0
0
b
0.2
.
0.8
0
0
0
c
0.8
0.2
.
0
0
0
d
0
0
0
.
0.6
0.4
e
0
0
0
0.4
.
0.6
f
0
0
0
0.6
0.4
.
;
u [*] :=
a 1
b 2
c 2
d 4
e 6
f 5
;
e) Come si pu`
o notare, rilassando l’integralit`a delle variabili x la soluzione ottima non `e pi`
u
intera.
La formulazione (M T Z) sembra apparentemente pi`
u semplice da risolvere in quanto pi`
u
compatta rispetto a (DF J) (presenta solo un numero polinomiale vincoli), tuttavia, come
confermato da questo esempio, il poliedro del rilassamento continuo di (M T Z) `e meno
stringente rispetto a quello del rilassamento continuo di (DF J).
Il confronto fra questi due poliedri viene analizzato in dettaglio nell’articolo The bad and the
good-and-ugly: formulations for the traveling salesman problem. di Gabor Pataki. In realt`
a
la formulazione (DF J) risulta essere l’unica che consente di trattare istanze di dimensione
elevata, perch`e, sebbene sia definita da un numero esponenziale di disuguaglianze, queste
ultime possono essere separate in modo estremamente efficiente, rendendo possibile un
Documento preparato da S. Coniglio, L. Taccari, V. Danese
13
laboratorio 2
Prof. E. Amaldi
Ottimizzazione
approccio con piani di taglio per risolvere il rilassamento continuo. Per risolvere il PLI all’ottimo si pu`
o adoperare un approccio Branch-and-Bound. Invece, la formulazione (M T Z)
ha un rilassamento continuo pi`
u debole rispetto a (DF J) che formisce bound nettamente
peggiori, rendendo pi`
u difficile l’applicazione del metodo di Branch-and-Bound.
La differenza tra le due formulazioni pu`o essere colta riscrivendo la formulazione (MTZ)
in funzione delle sole variabili xij . Utilizzando il Lemma di Farkas, si pu`o dimostrare che
la formulazione `e la seguente:
X
(MTZx ) min
cij xij
(16)
(i,j)∈E
s.t.
X
xij = 1
∀j ∈ V
(17)
xij = 1
∀i ∈ V
(18)
∀C ∈ C
(19)
∀ (i, j) ∈ E
(20)
i∈V
X
j∈V
X
xij ≤ |C| −
(i,j)∈C
|C|
|V | − 1
xij ∈ {0, 1}
dove C indica l’insieme di tutti i cicli orientati nel grafo (V \ {1} , E \ δ {1}) (identifichiamo
un ciclo C con una collezione di archi {(i, j)}).
I vincoli (19) sono implicati dai vincoli di eliminazione di sottociclo (SEC)
X
xij ≤ |S| − 1
∀S ⊂ V \ {1} , |S| ≥ 2
(i,j)∈δ(S)
Infatti, se, per esempio, consideriamo il sottoinsieme di nodi {2, 3, 4, 5} il vincolo (SEC) `e
x23 + x32 + x34 + x43 + x45 + x54 + x52 + x25 + x35 + x53 + x24 + x42 ≤ 3
(21)
Considerando la formulazione (MTZx ), i vincoli (19) scritti per tutti i sottocicli che passano
per i nodi {2, 3, 4, 5} sono
x23 + x34 + x45 + x52 ≤ 4 −
x32 + x25 + x54 + x43 ≤ 4 −
x23 + x35 + x54 + x42 ≤ 4 −
x32 + x45 + x53 + x24 ≤ 4 −
x34 + x25 + x53 + x42 ≤ 4 −
x43 + x52 + x35 + x24 ≤ 4 −
dove la quantit`
a 4−
4
|V |−1
4
|V | − 1
4
|V | − 1
4
|V | − 1
4
|V | − 1
4
|V | − 1
4
|V | − 1
(22)
(23)
(24)
(25)
(26)
(27)
`e sicuramente maggiore o uguale a 3 perch`e |V | = 6.
Documento preparato da S. Coniglio, L. Taccari, V. Danese
14

Download Report