laboratorio 6
Prof. E. Amaldi
Ottimizzazione
Metodo di Barriera (o dei Punti Interni) per la
Programmazione Lineare
Si consideri il seguente problema di programmazione lineare:
(P )
min
s.t.
cT x
(1)
Ax = b
(2)
x≥0
(3)
I vincoli che rendono difficile la soluzione del problema (P ) sono quelli di non negativit`
a sulle
variabili. Si consideri dunque l’applicazione del metodo di barriera logaritmica (metodo dei
punti interni) al problema in cui si rilassano solamente i vincoli x ≥ 0. Sia (P (µ)) il problema
che si ottiene per un dato valore del parametro di barriera µ, ovvero
(P (µ))
min
cT x − µ
n
X
ln xi
(4)
i=1
s.t.
Ax = b.
(5)
a) Si scrivano le condizioni di Karush-Kuhn-Tucker per il problema originale (P ). Le condizioni di KKT sono necessarie e/o sufficienti per l’ottimalit`a in ogni punto della regione
ammissibile?
Si considerino ora le condizioni di KKT del problema (P ): si perturbi ciascuna condizione di
complementariet`a con lo stesso parametro µ > 0. Si ottengono quindi le seguenti condizioni:
c − v + AT u = 0
(6)
Ax = b
(7)
x≥0
(8)
v≥0
(9)
xj v j = µ
per j = 1, . . . , n.
(10)
Definizione: Sia (x (µ) , u (µ) , v (µ)) una terna di vettori che soddisfi le condizioni (6)-(10).
Al variare del parametro µ > 0, le terne di vettori (x (µ) , u (µ) , v (µ)) descrivono una traiettoria detta cammino centrale (central path). Osserviamo che, se µ tende a 0, le condizioni di
complementariet`a si avvicinano a essere soddisfatte.
Pi`
u precisamente, si pu`o dimostrare che, se µ → 0, la terna (x (µ) , u (µ) , v (µ)) converge a
una terna (x∗ , u∗ , v ∗ ) che soddisfa le condizioni di KKT del problema (P ).
In particolare, si pu`o verificare che vale
L (x (µ) , u (µ) , v (µ)) = cT x (µ) − nµ ≤ cT x∗ ≤ cT x (µ)
dove x∗ `e una soluzione ottima per il problema (P ) (per la dimostrazione si veda l’Appendice 1).
Pertanto, lungo il cammino centrale, per ogni µ > 0, il valore della funzione obiettivo di (P )
u nµ dal valore ottimo cT x∗ .
cT x (µ) dista al pi`
Documento preparato da S. Coniglio, L. Taccari, V. Danese
1
laboratorio 6
Ottimizzazione
Prof. E. Amaldi
b) Si consideri ora il problema con barriera logaritmica (P (µ)) per un dato µ > 0. Si scrivano le condizioni di Karush-Kuhn-Tucker per tale problema. Le condizioni di KKT sono
necessarie e/o sufficienti per l’ottimalit`a in ogni punto della regione ammissibile?
c) Sia x∗ (µ) una soluzione che soddisfa le condizioni di KKT del problema (P (µ)) per determinati moltiplicatori u∗ (µ). La coppia (x∗ (µ) , u∗ (µ)) pu`o essere parte di una terna
(x∗ (µ) , u∗ (µ) , v ∗ (µ)) che appartenga al cammino centrale?
d) Si consideri l’algoritmo iterativo in cui si applica il metodo di Newton per la soluzione del
problema con barriera ad ogni iterazione. Si mostri come adattare la regola di derivazione
del passo di Newton in modo che le soluzioni soddisfino un insieme di vincoli lineari di
uguaglianza.
e) Si consideri il seguente problema:
min
x1 − x2
(11)
s.t. − x1 + x2 ≤ 1
(12)
x1 + x2 ≤ 3
(13)
x1 , x2 ≥ 0.
(14)
Si implementi il metodo in Matlab per risolverlo, partendo dal punto ammissibile x
¯ = (1, 1).
f) Si risolva il problema (11)-(14) con il metodo del simplesso, usando AMPL e il risolutore
CPLEX. Si verifichi che la soluzione trovata dal metodo dei punti interni non `e su un vertice
del poliedro ammissibile, al contrario di quella trovata da CPLEX. Qual’`e il motivo?
g) Si risolva ora il problema ottenuto modificando la funzione obiettivo in 21 x1 − x2 . Si
verifichi che la soluzione ottenuta col metodo dei punti interni coincide con quella trovata
dal metodo del simplesso. Per quale motivo?
Uno scheletro della funzione che implementa il metodo dei punti interni si trova in ipm stub.m:
Documento preparato da S. Coniglio, L. Taccari, V. Danese
2
laboratorio 6
Ottimizzazione
Prof. E. Amaldi
% interior point method for LPs of the form min c’*x : A*x <= b
function [xstar, vstar] = ipm(c, A, b, mu, x_init, gamma, epsilon)
OPTIONS = [];
[m, n] = size(A);
%%bring to standard form
s = b - A*x_init
x = [ x_init ; s ]
A = [A eye(m,m)]
c = [c zeros(1,m)]
%[m, n] = size(AA)
grad = zeros(m+n, 1);
H = zeros(m+n, m+n);
d = zeros(m+n, 1);
z = zeros(m, 1);
xks = [x_init’]
%% iterative method
while n * mu >= epsilon
%% compute gradient and Hessian
%% compute direction d with adapted Newton update
alpha = fminbnd(@(alpha) % write here the objective function f(x+alpha*d), ...
0, 1, OPTIONS);
xstar = x + alpha * d;
mu = gamma * mu;
x = xstar;
xks = [xks; x(1:n)’]
end
polyhedron_print(A,b); hold on;
plot(xks(:,1), xks(:,2), ’r.’);
end % end of function
Per rappresentare la regione ammissibile e la sequenza di soluzioni trovate, `e inclusa la
Documento preparato da S. Coniglio, L. Taccari, V. Danese
3
laboratorio 6
Ottimizzazione
Prof. E. Amaldi
funzione polyhedron print.m:
function polyhedron_print(A, b)
[m, n] = size(A);
Ineq = A(:, 1 : n-m);
rel = ’’;
for i = 1 : m
rel = [ rel ’<’ ];
end
t = extrpts(Ineq,rel,b);
t = t(1:2,:);
t = delcols(t);
t1 = t(1,:);
t2 = t(2,:);
z = convhull(t1,t2);
hold on
patch(t1(z),t2(z),[0.9, 0.9, 0.5])
end % end of function
che si appoggia agli script vr.m
function e = vr(m,i)
% The ith coordinate vector e in the m-dimensional Euclidean space.
e = zeros(m,1);
e(i) = 1;
delcols.m
function d = delcols(d)
% Delete duplicated columns of the matrix d.
d = union(d’,d’,’rows’)’;
n = size(d,2);
j = [];
for k =1:n
c = d(:,k);
for l=k+1:n
if norm(c - d(:,l),’inf’) <= 100*eps
j = [j l];
end
end
end
if ~isempty(j)
j = sort(j);
d(:,j) = [];
end
Documento preparato da S. Coniglio, L. Taccari, V. Danese
4
laboratorio 6
Ottimizzazione
Prof. E. Amaldi
e extrpts.m
function vert = extrpts(A, rel, b)
% Extreme points vert of the polyhedral set
%
X = {x: Ax <= b or Ax >= b, x >= 0}.
% Inequality signs are stored in the string rel, e.g.,
% rel = ’<<>’ stands for <= , <= , and >= , respectively.
[m, n] = size(A);
nlv = n;
for i=1:m
if(rel(i) == ’>’)
A = [A -vr(m,i)];
else
A = [A vr(m,i)];
end
if b(i) < 0
A(i,:) = - A(i,:);
b(i) = -b(i);
end
end
warning off
[m, n]= size(A);
b = b(:);
vert = [];
if (n >= m)
t = nchoosek(1:n,m);
nv = nchoosek(n,m);
for i=1:nv
y = zeros(n,1);
x = A(:,t(i,:))\b;
if all(x >= 0 & (x ~= inf & x ~= -inf))
y(t(i,:)) = x;
vert = [vert y];
end
end
else
error(’Number of equations is greater than the number of variables’)
end
vert = delcols(vert);
vert = vert(1:nlv,:);
Documento preparato da S. Coniglio, L. Taccari, V. Danese
5
laboratorio 6
Prof. E. Amaldi
Ottimizzazione
Appendice 1: calcolo del valore della lagrangiana di (P ) lungo il cammino
centrale
Sia (˜
x (µ) , u
˜ (µ) , v˜ (µ)) un punto appartenente al cammino centrale. In questo punto la
lagrangiana di (P ) vale:
L (˜
x (µ) , u
˜ (µ) , v˜ (µ)) = cT x
˜−
n
X
v˜j x
˜j + u
˜ T (A˜
x − b) = cT x
˜−
j=1
n
X
v˜j x
˜ j = cT x
˜ − nµ,
j=1
dato che (6)-(10) sono soddisfatte.
Inoltre, possiamo anche scrivere:
˜ (µ) , v˜ (µ)) = cT x
˜−
L (˜
x (µ) , u
n
X
v˜j x
˜j = −AT u
˜ + v˜
T
j=1
x−
n
X
T
˜ = −˜
uT A˜
uT b
v˜j x
˜j = (−AT u
˜) x
x = −˜
j=1
La lagrangiana di (P ) valutata in un punto del cammino centrale coincide quindi con il valore
della funzione obiettivo del problema duale.
[ Il duale di (P ) ha la seguente formulazione:
(D)
max bT u
s.t. AT u ≤ c.
Poich`e le variabili u sono libere da vincoli di non negativit`
a, il seguente problema ha lo stesso
valore ottimo del problema (D):
max − bT u
s.t. − AT u ≤ c
]
Indicando con x∗ una soluzione ottima per il problema (P ), valgono quindi le relazioni
cT x
˜ − nµ = −˜
uT b ≤ c T x ∗
T
∗
(per la dualit`a debole in programmazione lineare)
(perch´e cT x∗ `e il valore ottimo),
T
˜
c x ≤c x
che implicano
cT x
˜ − nµ ≤ cT x∗ ≤ cT x
˜
Pertanto, lungo il cammino centrale, facendo convergere µ a 0, la lagrangiana di (P ) tende al
valore ottimo della funzione obiettivo.
Documento preparato da S. Coniglio, L. Taccari, V. Danese
6
laboratorio 6
Prof. E. Amaldi
Ottimizzazione
Soluzione
a) Siano u e v, rispettivamente, i moltiplicatori dei vincoli di uguaglianza e di non negativit`
a.
La funzione Lagrangiana di P `e:
L(x, v, u) = cT x −
n
X
vj xj + uT (Ax − b)
j=1
e le condizioni KKT sono:
c − v + AT u = 0
(15)
Ax = b
(16)
x≥0
(17)
v≥0
(18)
xj v j = 0
per j = 1, . . . , n.
(19)
Si ricordi che la condizione di qualifica dei vincoli `e sempre soddisfatta nel caso in cui
i vincoli siano lineari, anche nel caso in cui questi non abbiano gradienti linearmente
indipendenti in ogni punto della regione ammissibile (si faccia riferimento, ad esempio, al
libro Numerical Optimization di J. Nocedal e S. Wright, pagina 364). Quindi le KKT sono
necessarie.
Inoltre, si osservi che il problema (P ) `e convesso, in quanto sia la funzione obiettivo, sia
i vincoli, sono lineari. Quindi la funzione lagrangiana di (P ) `e convessa e le condizioni di
KKT sono anche sufficienti.
b) La funzione Lagrangiana di P (µ) `e:
Lµ (x, u) = cT x − µ
n
X
ln xj + uT (Ax − b)
j=1
e le condizioni KKT, ponendo β = β(x) :=
1
1
,...,
x1
xn
!T
, sono
c − µβ + AT u = 0
(20)
Ax = b.
(21)
Si osservi che non compare esplicitamente la condizione x ≥ 0 ma ogni soluzione x di P (µ)
in cui la funzione obiettivo assume valore finito deve soddisfare x > 0. La condizione di
qualifica dei vincoli `e sempre soddisfatta perch`e i vincoli sono lineari (si faccia riferimento
al punto a)). Quindi le KKT sono necessarie.
Inoltre, si osservi che il problema (P (µ)) `e convesso, in quanto la funzione obiettivo `e data
dalla somma di un termine lineare con una funzione convessa moltiplicata per un parametro
positivo e i vincoli sono lineari. Quindi la funzione lagrangiana di (P ) `e convessa e le KKT
sono anche sufficienti.
c) Sia x∗ , u∗ una soluzione di (20)-(21).
dev’essere:
Affinch`e la coppia di vettori soddisfi (6)-(10),
!T
1
1
v ∗ (µ) = µ
,...,
,
(22)
x1
xn
Documento preparato da S. Coniglio, L. Taccari, V. Danese
7
laboratorio 6
Prof. E. Amaldi
Ottimizzazione
da (22) deriviamo:
xj vj = µ.
(23)
La terna (x∗ , u∗ , v ∗ ) cos`ı definita appartiene dunque al cammino centrale. Pertanto, quando µ → 0, la soluzione ottima del problema (P (µ)) converge a una soluzione ottima del
problema (P ).
¯ `e data da:
d) Per il problema min f (x), la direzione di Newton in un punto x
d = −(∇2 f (¯
x))−1 ∇f (¯
x).
Per un problema con vincoli di uguaglianza del tipo:
min f (x) : Ax = b,
(24)
dobbiamo trovare una direzione di discesa d tale che il nuovo punto x + d, partendo da
un x che soddisfi Ax = b, sia anch’esso ammissibile, soddisfando A(x + d) = b. Questo
implica Ad = 0.
Le condizioni KKT di (24) sono:
∇f (x) + AT u = 0
(25)
Ax = b.
(26)
Applicandole al problema
min f (x + d) : A(x + d) = b
(27)
∇f (x + d) + AT u = 0
(28)
Ad = 0.
(29)
d
abbiamo:
Applichiamo ora il metodo di Newton, approssimando in serie di Taylor al primo ordine il
sistema nonlineare nell’intorno di d = 0. Otteniamo:
∇f (x) + ∇2 f (x)d + AT u = 0
(30)
Ad = 0
(31)
che, riscritto in forma matriciale, diventa:
!
∇2 f (x) AT
A
0
d
u
!
=
−∇f (x)
0
!
.
L’aggiornamento iterativo di x `e xk+1 = xk + αk dk , dove αk `e ottenuto mediante ricerca
unidimensionale. Se tale ricerca `e esatta, abbiamo:
αk = argminα≥0 f (xk + αdk ).
Si noti che, se xk `e ammissibile, anche xk + αk dk lo `e. Infatti, Axk+1 = Axk + Aαk dk = b,
dato che Adk = 0 per costruzione.
Documento preparato da S. Coniglio, L. Taccari, V. Danese
8
laboratorio 6
Prof. E. Amaldi
Ottimizzazione
e) Portiamo il problema originale in forma standard:
(P )
min
x1 − x2
s.t.
− x1 + x2 + x3 = 1
x1 + x2 + x4 = 3
x1 , x2 , x3 , x4 ≥ 0.
Il corrispondente problema con barriera logaritmica `e:
(P (µ))
min
x1 − x2 − µ
4
X
ln xj
j=1
s.t.
− x1 + x2 + x3 = 1
x1 + x2 + x4 = 3.
La matrice A `e:
−1 1 1 0
1
1 0 1
!
.
Il gradiente della funzione obiettivo `e:
µ
µ
µ
µ T
1 − , −1 − , − , −
x1
x2 x3 x4
e l’hessiano `e:
diag
µ µ µ µ
, , ,
x21 x22 x23 x24
.
Il sistema da risolvere nel metodo di Newton `e quindi:
µ
x21
0
0
0
−1 1
0
µ
x22
0
0
1
0
0
µ
x23
0
1
0
0
0
0
µ
x24
−1
1
1
0
0
1
1
0
1
0
1
0
1
0
0
d1
µ
x1
µ
x2
−1 +
d2
1+
µ
d3
x3
=
µ
d4
x4
u1
0
u2
0
,
Sia ǫ > 0, γ < 1, k = 0, µ0 > 0 e x(µ0 ) ammissibile per (P ), ossia tale che x(µ0 ) > 0 e
Ax(µ0 ) = b. Lo schema iterativo del metodo `e il seguente:
1. Si risolve P (µk ) partendo dal punto x(µk ), ottenendo una soluzione x(µk )∗ .
2. Se nµk < ǫ, l’algoritmo termina con soluzione x∗ (µk ).
3. Si aggiornano µk+1 := γµk , x(µk+1 ) := x(µk )∗ e k := k + 1.
4. Si ripete dal passo 1.
Documento preparato da S. Coniglio, L. Taccari, V. Danese
9
laboratorio 6
Ottimizzazione
Prof. E. Amaldi
Il file Matlab che implementa il metodo `e ipm.m
% interior point method for LPs of the form min c’*x : A*x <= b
function [xstar, vstar] = ipm(c, A, b, mu, x_init, gamma, epsilon)
OPTIONS = [];
[m, n] = size(A);
%%bring to standard form
s = b - A*x_init
x = [ x_init ; s ]
A = [A eye(m,m)]
c = [c zeros(1,m)]
%[m, n] = size(AA)
grad = zeros(m+n, 1);
H = zeros(m+n, m+n);
d = zeros(m+n, 1);
z = zeros(m, 1);
xks = [x_init’]
%% iterative method
while n * mu >= epsilon
for i = 1 : m+n
grad(i) = c(i) - mu / x(i)
H(i,i) = mu/x(i)^2;
end
H
A
N = [ H, A’;
A, zeros(m, m) ]
bN = [ -grad; zeros(m, 1) ]
xN = N \ bN
d = xN(1 : m+n)
z = xN(m+n + 1 : m+n + m)
alpha = fminbnd(@(alpha) c*(x + alpha*d) - mu*sum(log(x + alpha*d),1), ...
0, 1, OPTIONS);
xstar = x + alpha * d;
mu = gamma * mu;
x = xstar;
xks = [xks; x(1:n)’]
end
vstar = mu ./ xstar;
polyhedron_print(A,b); hold on;
plot(xks(:,1), xks(:,2), ’r.’);
end % end of function
Utilizziamo i parametri µ0 = 1, γ = 0.9, ǫ = 0.001. Si ricordi che c = (1, −1, 0, 0)T ,
b = (1, 3)T e x(µ0 ) = (1, 1)T .
Documento preparato da S. Coniglio, L. Taccari, V. Danese
10
laboratorio 6
Prof. E. Amaldi
Ottimizzazione
Per eseguire il codice usiamo il comando:
ipm([1 -1 ], [-1 1
; 1 1], [1; 3], 1, [1; 1], 0.5, 0.01)
ottenendo la soluzione x∗ = (0.58, 1.58, 0.00, 0.85)T , v = (0.00, 0.00, 0.50, 0.00)T . L’algoritmo traccia la successione di punti illustrata in Figura 1. Si noti che la soluzione trovata
non `e su un vertice del poliedro.
2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0
0
0.5
1
1.5
2
2.5
3
Figura 1: Punti x(µk ) individuati dall’algoritmo.
f) Utilizziamo il modello AMPL ipm-simplex.mod
# file: ipm-simplex.mod
var x{1..4} >= 0;
minimize f : 0.5*x[1] - x[2];
c1 : -x[1] + x[2] + x[3] = 1;
c2 : x[1] + x[2] + x[4] = 3;
e il file ipm-simplex.run
# file: ipm-simplex.run
model ipm-simplex.mod
option solver cplex;
solve;
display x;
display x.rc;
Risolvendo con CPLEX, otteniamo la soluzione xs = (0, 1, 0, 2)T .
L’istruzione display x.rc stampa i costi ridotti delle variabili nella soluzione ottima
trovata. Il fatto che vi siano costi ridotti nulli indica che la soluzione ottima non sia
unica (`e possibile effettuare un cambio di base, facendo pivot su una delle variabili a costo
ridotto nullo, ottenendo un’altra soluzione di costo uguale). Analiticamente, questo `e
conseguenza del fatto che il gradiente della funzione obiettivo `e perpendicolare alla retta
−x1 + x2 = 1, corrispondente al primo vincolo. Il poliedro ha quindi una faccetta ottima,
quella congiungente i punti del segmento tra R = (0, 1) e S = (1, 2). Tutti i punti nel
segmento hanno lo stesso valore di funzione obiettivo, ma un diverso valore del termine di
barriera. La soluzione individuata dal metodo dei punti interni, detta centro analitico, `e
una soluzione che, tra quelle ottime, minimizza anche il termine di barriera.
Documento preparato da S. Coniglio, L. Taccari, V. Danese
11
laboratorio 6
Prof. E. Amaldi
Ottimizzazione
T
g) Risolviamo con i due metodi il problema modificato in cui c = ( 12 , −1, 0, 0) , ottenendo
x∗ = (1.000, 2.000, 0.000, 0.000)
v ∗ = (0.000, 0.000, 0.6750, 0.2250)
xs = (1, 2, 0, 0)
v s = (0, 0, 0.75, 0.25),
` evidente che la soluzione del metodo dei punti interni
con un tolleranza ǫ = 10−5 . E
converge allo stesso valore di quella del simplesso, come si pu`o osservare dal cammino
centrale (central path)in Figura 2. Dato che il gradiente della funzione obiettivo non `e
parallelo ad alcun vincolo, il problema ammette una singola soluzione ottima.
2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0
0
0.5
1
1.5
2
2.5
3
Figura 2: Punti x(µk ) individuati dall’algoritmo applicato al problema modificato.
(Facoltativo) Appendice 2: confronto con il metodo del simplesso
Consideriamo l’algoritmo del simplesso primale applicato al problema (P ). Riprendiamo le KKT
del problema (P ) e vediamo quali vengono violate prima che tale algoritmo termini:
c − v + AT u = 0
(32)
Ax = b
(33)
x≥0
(34)
v≥0
(35)
xj v j = 0
per j = 1, . . . , n.
(36)
Il metodo del simplesso propone a ogni iterazione k una soluzione ammissibile per il primale (P ),
ossia tale da soddisfare Axk = b e xk ≥ 0. Sia Bk la sottomatrice di A formata dalle colonne
associate alle variabili in base all’iterazione k, Nk la sottomatrice corrispondente alle variabili
fuori base. Abbiamo dunque
A = [Bk | Nk ]
e
xk = (xBk , xNk ) = (Bk −1 b, 0)
Documento preparato da S. Coniglio, L. Taccari, V. Danese
12
laboratorio 6
Prof. E. Amaldi
Ottimizzazione
Ad ogni iterazione k, `e possibile associare alla soluzione corrente xk un vettore (uk , v k ) costruito
in modo tale da soddisfare le condizioni (32) e (36). Sia xkj la j-esima componente del vettore
xk . Per garantire il rispetto della complementariet`a, innanzitutto imponiamo
vkj = 0 se xkj `e in base.
Per garantire la condizione (32), imponiamo
AT uk + v k = ck ,
ovvero
Bk T u k + v B k = c B k
N k T u k + v Nk = cNk .
Poich`e abbiamo imposto v Bk = 0, otteniamo
uk = (Bk T )
−1
c Bk
v Nk = cNk − Nk T (Bk T )
−1
cB k
L’unica condizione di KKT che viene violata prima che l’algoritmo termini `e quindi
v Nk ≥ 0.
In particolare, osserviamo che il vettore v Nk coincide esattamente con i costi ridotti delle variabili
fuori base: l’algoritmo termina quando tutti i costi ridotti sono non negativi, ovvero quando
v Nk ≥ 0.
Documento preparato da S. Coniglio, L. Taccari, V. Danese
13
© Copyright 2024 ExpyDoc
