Complementi di Ricerca Operativa 6.1 Controllo ottimo dell

Esercitazione 6
Complementi di Ricerca Operativa
Prof. E. Amaldi
6.1 Controllo ottimo dell’instradamento e del traffico in rete. Si consideri un
grafo diretto G = (V, A), che modella una rete di telecomunicazione, ed una matrice
D di dimensione |V | × |V |, in cui ogni elemento (dij ) rappresenta la quantità di dati
che deve essere inviata dal nodo i al nodo j (i e j rappresentano rispettivamente
l’origine e la destinazione del flusso di traffico). L’instradamento di una quantità
di dati f lungo ogni arco (i, j) comporta un costo cij (f ) (ritardo di propagazione,
probabilità di perdita di dati ecc.).
(a) Proporre un modello di PNL per individuare l’instradamento ottimo del traffico
sulla rete (suggerimento: si utilizzi un insieme di variabili xp , ad ognuna delle
quali corrisponde un cammino sul grafo G).
(b) Assumendo che le funzioni di costo cij (f ) siano convesse e continuamente
differenziabili, che relazione lega i valori delle derivate direzionali
∂cij (f )/∂xp
ai valori assunti dalle variabili in una soluzione ottima?
6.2 Metodi di barriera in PL. Si consideri il problema di PL in forma standard:
[P]
min{c⊤ x | Ax = b, x ≥ 0}.
Sia x∗ la soluzione di P. Si consideri ora il seguente problema, in cui è stata applicata
una funzione di barriera ai vincoli di non negatività:
[P(µ)]
min c⊤ x − µ
n
P
ln xj
j=1
s.t. Ax = b.
Sia x∗ (µ) la soluzione di P(µ). Si verifichi che x∗ (µ) converge a x∗ per µ → 0.
[Suggerimento: si confrontino le condizioni di ottimalità di P e P(µ).]
6.3 Metodo di penalità quadratica e dei Lagrangeani aumentati.
(a) Considerare il problema di programmazione non lineare:
min
x21 + x22
s.v. x1 + x2 = 1
e determinarne la soluzione ottima con il rispettivo moltiplicatore di Lagrange.
Come varia la soluzione ottima del problema di penalità quadratica quando il
parametro di penalità µ tende a 0? Calcolare la matrice Hessiana del problema
di penalità e i suoi autovalori in funzione di µ. Tracciare le linee di livello
della funzione di penalità quadratica Q(x, µ) per µ = 1 e µ = 0.1. Cosa può
succedere se il problema di penalità è risolto mediante il metodo del gradiente?
Eseguire tre iterazioni del metodo basato sulla Lagrangiana aumentata (metodo
dei moltiplicatori) partendo da µ = 1 e v = 0. Cosa si osserva?
Documento preparato da E. Amaldi e A. Ceselli
1
Esercitazione 6
Complementi di Ricerca Operativa
Prof. E. Amaldi
(b) Risolvere il problema:
min (x1 − 2)4 + (x1 − 2x2 )2
s.v.
x21 − x2 = 0
mediante il metodo di penalità quadratica, partendo dal punto iniziale x0 =
(2, 1) e con µ = 10..
Rappresentare nel piano x1 , x2 le linee di livello della funzione obiettivo, la
curva relativa al vincolo e i primi 5 punti generati dall’algoritmo.
6.4 Metodi di penalità e dei Lagrangiani aumentati. Si consideri il problema di
programmazione non lineare:
min
x1 x2
s.v. x1 − 2x2 = 3
(a) Determinare analiticamente la soluzione ottima.
(b) Per quale valore del parametro di penalità la funzione di penalità quadratica
ha un minimo? Determinare il minimo in funzione del parametro di penalità
µ. Qual’è il limite?
(c) Per un valore sufficientemente elevato del parametro di penalità, applicare il
metodo basato sulla Lagrangiana aumentata, partendo dal valore del moltiplicatore di Lagrange v = 0.
(d) Formulare il problema di penalità esatta e determinare il minimo non vincolato,
se esiste.
6.5 Metodo di programmazione quadratica sequenziale. Si consideri il problema
di programmazione non lineare:
min
x1 − x2
2
s.v. x1 + x22 ≤ 1
(a) Partendo dal punto (−1, 0) e dal valore iniziale del moltiplicatore di Lagrange
u = 1, eseguire due iterazioni del metodo di programmazione quadratica sequenziale. Spiegare perché viene scelto alla prima iterazione un valore positivo
del moltiplicatore.
6.6 Risolvere con Matlab il problema di ottimizzazione non lineare:
min exp(x1 x2 x3 x4 x5 ) − 12 (x31 + x32 + 1)2
s.v.
x21 + x22 + x23 + x24 + x25 = 10
x2 x3 − 5x4 x5 = 0
x31 + x32 = −1
partendo dal punto iniziale x0 = (−1.8, 1.7, 1.9, −0.6, −0.8).
Documento preparato da E. Amaldi e A. Ceselli
2
Esercitazione 6
Complementi di Ricerca Operativa
Prof. E. Amaldi
Soluzioni
6.1 Controllo ottimo dell’instradamento e del traffico in rete.
(a) Sia Pij l’insieme dei cammini sul grafo G(V, A) che partono da i ∈ V (sorgente)
ed arrivano in j ∈ V (destinazione). Tutti gli archi in questi cammini sono
orientati nella direzione dalla sorgente alla destinazione. Per ogni cammino
p ∈ Pij , sia xp ∈ R+ una variabile che rappresenta la porzione del traffico dij
tra i e j instradata lungo il cammino p. Per ogni coppia origine–destinazione
(i, j), l’intero flusso di traffico tra i e j deve essere instradato:
X
xp = dij ∀(i, j) ∈ V × V.
p∈Pij
Il flusso totale fij su ogni arco (i, j) ∈ A può essere calcolato come:
X
fij =
xp
cammini p che utilizzano (i,j)
Esistono due importanti tipi di instradamento ottimo: l’ottimo per l’utente, in
cui ogni utente cerca indipendentemente la strada per lui più breve, e l’ottimo
per il sistema, in cui l’obiettivo è minimizzare il tempo medio di percorrenza su
tutti gli utenti, eventualmente sfavorendone alcuni. Cerchiamo di determinare
l’ottimo per il sistema: il costo di ogni soluzione è
X
cij (fij )
D(x) =
(i,j)∈A
Il problema può essere, quindi, formulato come segue:
X
X
minimize D(x) =
cij (
xp )
(i,j)∈A
s.t.
X
cammini p che utilizzano (i,j)
xp = dij
∀(i, j) ∈ V × V
p∈Pij
xp ≥ 0
∀(i, j) ∈ V × V, ∀p ∈ Pij .
Questo problema è noto come flusso multiprodotto. In Figura 1 è riportato un
esempio di soluzione: il traffico tra il vertice 1 ed il vertice 5 viene instradato
lungo 3 cammini:
A= 1→2→5
(2 Mb)
B = 1 → 2 → 4 → 5 (1 Mb)
C = 1 → 3 → 4 → 5 (1.5 Mb)
quindi xA = 2, xB = 1 e xC = 1.5. E’ possibile individuare altri cammini
tra il vertice 1 ed il vertice 5, ad esempio: D = 1 → 2 → 3 → 2 → 5. La
corrispondente variabile xD assume il valore 0. Il traffico instradato lungo gli
archi è riportato nella seguente tabella:
Documento preparato da E. Amaldi e A. Ceselli
3
Complementi di Ricerca Operativa
Esercitazione 6
Prof. E. Amaldi
Figura 1: Esempio di instradamento di traffico tra due vertici nell’esercizio 6.1.
(1,2)
(2,5)
(2,4)
(1,3)
(3,4)
(4,5)
3Mb
2Mb
1Mb
1.5Mb
1.5Mb
2.5Mb
(b) Le condizioni di ottimalità del primo ordine affermano che un punto x∗ è il
minimo (globale) di una funzione strettamente convessa f su X se e solo se
∇f (x∗ )(x − x∗ ) ≥ 0 ∀x ∈ X.
Innanzitutto, riportiamo un’osservazione inerente i problemi di ottimizzazione
non lineare vincolata nei quali l’insieme delle soluzioni ammissibili X è un
simplesso:
n
X
X = {x|x ≥ 0,
xi = r}.
i=1
In questo caso, secondo le condizioni del primo ordine, x∗ = (x∗1 , . . . , x∗n ) è un
punto di minimo solo se
n
X
∂f (x∗ )
i=1
∂xi
(xi −
x∗i )
≥ 0, ∀xi ≥ 0,
n
X
xi = r
i=1
Ora, si fissi un indice i per il quale x∗i > 0 e sia j un qualsiasi altro indice.
Considerando un punto ammissibile x = (x1 , . . . , xn ) con xi = 0, xj = x∗j + x∗i ,
xm = x∗m ∀m 6= i, j, e sostituendo nella precedente espressione, si ottiene
∂f (x∗ ) ∂f (x∗ )
x∗i ≥ 0
−
∂xj
∂xi
o, in modo equivalente,
x∗i > 0 implica che
∂f (x∗ )
∂f (x∗ )
≤
∂xi
∂xj
∀j.
(1)
E’ possibile dimostrare che, nel caso la funzione f sia convessa, questa condizione è anche sufficiente per l’ottimalità [Bert03]; inoltre questo risultato
Documento preparato da E. Amaldi e A. Ceselli
4
Esercitazione 6
Complementi di Ricerca Operativa
Prof. E. Amaldi
può essere generalizzato al caso in cui l’insieme delle soluzioni ammissibili non
sia un simplesso, ma il prodotto cartesiano di più simplessi (una condizione
analoga vale per ognuno dei simplessi nel prodotto cartesiano [Bert03]).
L’insieme delle soluzioni ammissibili per il problema di flusso in esame è esattamente un prodotto cartesiano di simplessi (uno per ogni coppia originedestinazione). La derivata prima della funzione costo rispetto alle variabili
xp è dunque:
X
∂cij (f )
∂D(x)
=
|
∂xp
∂f fij
archi (i,j) sul cammino p
ed, utilizzando la condizione (1), possiamo affermare che un insieme di cammini
di flusso è ottimo se e solo se il flusso è instradato solo lungo cammini p aventi
∂D(x)
minima.
∂xp
Per cogliere l’idea più intuitivamente, supponiamo che tij (f ) rappresenti la lunghezza (o tempo di percorrenza) dell’arco (i, j) quando lungo l’arco è instradata
una quantità di flusso f . Poniamo
Z fij
tij (f ) df
cij (fij ) =
0
Ora, in ogni soluzione ottima per il sistema rispetto alle funzioni cij (f ) il flusso
tra ogni coppia origine-destinazione (i, j) è instradato esclusivamente lungo
cammini da i a j a lunghezza minima (e aventi la stessa lunghezza). Ovvero, a
patto di definire le funzioni cij (f ) in modo opportuno, è possibile individuare
l’ottimo per l’utente cercando l’ottimo per il sistema.
6.2 Metodi di barriera in PL. Mostriamo prima di tutto che per ogni soluzione x∗ (µ)
(ottima in P(µ) e ammissibile in P) si può facilmente calcolare la corrispondente
soluzione duale ammissibile (sia in P(µ) che in P). Si considerino le variabili duali
uj relative ai vincoli di nonnegatività di x espressi nella forma “inversa” −xj ≤ 0
per j ≤ n. La funzione Lagrangiana del problema P è:
⊤
L(x, u, v) = c x −
n
X
uj xj + v(Ax − b),
j=1
e la funzione Lagrangiana del problema P(µ) è:
⊤
Lµ (x, v) = c x(µ) − µ
n
X
ln xj (µ) + v(Ax − b).
j=1
Richiedere che x sia primale ammissibile e (u, v) duale ammissibile di P equivale a
imporre la condizione1 ∇L = 0, ovvero
∀j ≤ n (cj − uj + vAj = 0),
1
Nella teoria dell’algoritmo del simplesso, l’ammissibilità primale e duale implica l’ottimalità. Questo
succede perché le soluzioni di base sono sempre ai vertici del poliedro; considerando i punti interni al
poliedro, si può avere l’ammissibilità primale e duale senza l’ottimalità.
Documento preparato da E. Amaldi e A. Ceselli
5
Esercitazione 6
Complementi di Ricerca Operativa
Prof. E. Amaldi
dove Aj è la j-esima colonna di A. Ora imponiamo che x sia primale ammissibile e
v duale ammissibile in P(µ) mediante l’equazione ∇Lµ = 0, ovvero
∀j ≤ n (cj −
µ
+ vAj = 0).
xj
Si consideri ora che essendo x∗ (µ) una soluzione di P(µ), soddisfa la condizione di
ottimalità ∇Lµ (x∗ (µ), v ∗ ) = 0 (ammissibilità primale/duale in P(µ)). Definendo
u∗j (µ) = x∗µ(µ) otteniamo ∇L(x∗ (µ), u∗ (µ), v ∗ ) = 0, ovvero l’ammissibilità primaj
le/duale in P, perciò x∗ (µ) è primale ammissibile e u∗ (µ) duale ammissibile (corrispondente ai vincoli x ≥ 0) del problema P. In sintesi, da ogni soluzione ottima
x∗ (µ) di P(µ) possiamo facilmente ricavare una soluzione duale ammissibile u∗ (µ)
di P tale che x∗j (µ)u∗j (µ) = µ per ogni j ≤ n.
Ora, per µ → 0, otteniamo che la coppia primale/duale (x∗ (µ), u∗ (µ)) converge
a una coppia (x′ , u′ ). Si noti che dato che le condizioni cj − yj + zAj = 0 sono
lineari, che y ∗ (µ) le soddisfa e che y ∗ (µ) → y ′ , anche y ′ le soddisfa. Dato inoltre
che x∗j (µ)u∗j (µ) = µ per ogni j ≤ n e per ogni µ > 0, si ha che x∗j (µ)u∗j (µ) → 0 per
µ → 0 e quindi x′j u′j = 0 per ogni j ≤ n. Quindi per costruzione x′ e u′ soddisfano
le equazioni degli scarti complementari ∀j ≤ n (xj uj = 0) derivate dal problema P.
Ciò significa che (x′ , u′ ) è una coppia primale/duale ottima per P, ovvero x′ = x∗ ,
come volevasi dimostrare.
6.3 Metodo di penalità quadratica e dei Lagrangeani aumentati.
(a) La funzione Lagrangeana in questo caso è
L(x, v) = x21 + x22 + v(x1 + x2 − 1)
e, secondo le condizioni KKT, in ogni soluzione ottima
1 +v
=0
∇x L(x, v) = 2x
2x2 +v
x1 + x2 = 1
da cui si deriva facilmente che la soluzione ottima è x1 = x2 = 0.5 con v = −1.
La funzione con penalità quadratica Q(x, µ) è definita come segue:
Q(x, µ) = x21 + x22 +
1
(x1 + x2 − 1)2
2µ
e
∂Q(x, µ)
1
1
= (2 + )x1 + (x2 − 1)
∂x1
µ
µ
∂Q(x, µ)
1
1
= (2 + )x2 + (x1 − 1).
∂x2
µ
µ
La matrice Hessiana è quindi
H=
Documento preparato da E. Amaldi e A. Ceselli
2+
1
µ
1
µ
1
µ
2+
1
µ
6
Complementi di Ricerca Operativa
Esercitazione 6
0.8
0.8
0.6
0.6
x2
1
x2
1
Prof. E. Amaldi
0.4
0.4
0.2
0.2
0
0
0
0.2
0.4
0.6
0.8
1
0
0.2
x1
0.4
0.6
0.8
1
x1
Figura 2: funzione di penalità quadratica Q(x, µ) nell’esercizio 6.3 (a sinistra per µ = 1.0
ed a destra per µ = 0.1)
avente autovalori λ1 = 2 e λ2 = 2/µ + 1. Per µ che tende a 0 la differenza tra
gli autovalori tende all’infinito. Questo è un indicatore, ad esempio, di cattive
prestazioni per il metodo del gradiente quando µ è molto piccolo. Le curve di
livello di Q(x, µ) per µ = 1 e µ = 0.1 sono rappresentate in Figura 2.
La funzione Lagrangeana aumentata è definita come:
LA (x, v, µ) = x21 + x22 + v(x1 + x2 − 1) +
1
(x1 + x2 − 1)2
2µ
Nel metodo della Lagrangeana aumentata, ad ogni iterazione k si considerano µ e
v dei parametri fissati rispettivamente ai valori µk e v k , e si risolve (solitamente in
modo approssimato) il rimanente problema nelle variabili x, ottenendo una soluzione
(xk1 , xk2 ). Poi si aggiorna il valore del moltiplicatore v in funzione della violazione
nel corrispondente vincolo:
v k+1 = v k +
1
h(xk1 , xk2 )
k
µ
dove h(xk1 , xk2 ) = xk1 + xk2 − 1 è l’entità della violazione del vincolo (notare l’analogia
con i metodi del sottogradiente). Infine, si decrementa il valore del parametro µ
(solitamente, in dipendenza del comportamento durante l’ottimizzazione rispetto a
x). In questo esempio, per semplicità, dimezziamo il valore del parametro µ ad ogni
Documento preparato da E. Amaldi e A. Ceselli
7
Esercitazione 6
Complementi di Ricerca Operativa
Prof. E. Amaldi
iterazione. Nell’esempio in esame, si ottiene il comportamento seguente:
iter.
µ
v
x
1
1.0
0.0
(0.25, 0.25)
2
0.5
0.5
(0.25, 0.25)
3
0.25
1.0
(0.3, 0.3)
4
0.125 1.6 (0.3556, 0.3556)
5 0.0625 2.32 (0.4024, 0.4024)
(b) In questo caso, la funzione di penalità quadratica è:
Q(x, µ) = (x1 − 2)4 + (x1 − 2x2 )2 +
1 2
(x − x2 )2
2µ 1
As ogni iterazione k del metodo di penalità quadratica, si fissa a µk il valore del
parametro µ e si risolve (solitamente in modo approssimato) il rimanente problema
non vincolato nelle variabili x, ottenendo una soluzione xk . Poi si decrementa il
valore del parametro µ e si itera, utilizzando xk come nuovo punto iniziale. In
questo esempio, per semplicità, dividiamo per 10 il valore del parametro µ ad ogni
iterazione, ottenendo il seguente comportamento.
iter.
0
1
2
3
4
5
µ
x
(2, 1)
5.0
(1.4539, 0.7608)
0.5
(1.1687, 0.7407)
0.05
(0.9906, 0.8425) 0.005
(0.9508, 0.8875) 0.0005
(0.9456, 0.8941) 0.00005
le curve di livello della funzione di penalità ed i punti di minimo per le prime
iterazioni sono riportate in figura 3. A patto di scegliere un’adeguata politica per
l’aggiornamento dei moltiplicatori, il metodo converge molto rapidamente.
6.4 Metodi di penalità e dei Lagrangeani aumentati.
(a) Ponendo x1 = 3 + 2x2 , l’obiettivo diventa
min(3 + 2x2 ) · x2
e nell’unico punto di minimo il valore della derivata rispetto a x2 è 0:
3 + 4x2 = 0 → x2 = −3/4 e quindi x1 = 3/2
(b) La funzione di penalità quadratica è la seguente:
Q(x, µ) = x1 x2 +
Documento preparato da E. Amaldi e A. Ceselli
1
(x1 − 2x2 − 3)2 .
2µ
8
Esercitazione 6
Complementi di Ricerca Operativa
Prof. E. Amaldi
Figura 3: funzione di penalità quadratica Q(x, µ) nell’esercizio 6.3(b) durante le prime
iterazioni
Documento preparato da E. Amaldi e A. Ceselli
9
Esercitazione 6
Complementi di Ricerca Operativa
Prof. E. Amaldi
Studiandone le derivate parziali
1
∂Q
2(x1 − 2x2 − 3)
= x2 +
∂x1
2µ
e
1
∂Q
(−4)(x1 − 2x2 − 3)
= x1 +
∂x2
2µ
si può trovare il punto di minimo in funzione del parametro µ:
x1 − 2x2 − 3 = µ(−x2 )
x1 − 2x2 − 3 = µx2 1
quindi
x2 +
1
(−4x2
µ
x1 = −2x2
− 3) = 0 → (1 − µ4 )x2 =
3
µ
3
6
, µ−4
). La matrice
e la soluzione ottima in funzione del parametro µ è x∗µ (− µ−4
Hessiana è la seguente:
1
1 − µ1
µ
Hµ =
4
− µ2
µ
che è definita positiva per ogni µ > 0. Un minimo esiste per µ 6= 4. Per µ
tendente a 0:
3 3
lim x∗µ = ( , − )
µ→0
2 4
che è in effetti la soluzione ottima del problema.
(c) Lasciata agli studenti (analogo agli esercizi precedenti).
(d) Come funzione di penalità esatta utilizziamo la funzione valore assoluto:
1
|x1 − 2x2 − 3|
2µ
min x1 x2 +
Se x1 − 2x2 − 3 ≥ 0
x2 +
1
= 0 → x2 = − 2µ
1
2µ
x1 − 2 ·
1
2µ
= 0 → x1 =
1
µ
altrimenti,
1
x2 + − 2µ
= 0 → x2 =
x1 + 2 ·
1
2µ
1
2µ
= 0 → x1 = − µ1
Quindi,
(
1
( µ1 , − 2µ
)
x∗µ =
1
(− µ1 , 2µ
)
Documento preparato da E. Amaldi e A. Ceselli
se x1 − 2x2 − 3 ≥ 0
altrimenti.
10
Esercitazione 6
Complementi di Ricerca Operativa
Prof. E. Amaldi
Il termine di penalità è minimizzato nel punto di non differenziabilità del valore
assoluto:
1
2
1
+2·
−3≥0 → µ≤ .
µ
2µ
3
L’ottimo del problema è individuato per µ = 2/3.
6.5 Metodo di programmazione quadratica sequenziale.
L’idea alla base del metodo di programmazione quadratica sequenziale è la seguente:
come nei metodi delle tangenti (tipo Newton) per ottimizzazione non vincolata,
invece di risolvere il sistema delle condizioni KKT per il problema di partenza (che
può essere non lineare), ne costruisco un’approssimazione locale (lineare). Nel caso
in esame, scelto un punto x̄:
f (x) = x1 − x2 →
(∇f (x̄))T (x − x̄) = (1, −1)T (x1 − x¯1 , x2 − x¯2 )
e
g(x) = x21 − x22 − 1 →
(g(x̄) + ∇g(x̄))T (x − x̄) = (x¯1 2 + x¯2 2 − 1) + (2x1 , −2x2 )T (x1 − x¯1 , x2 − x¯2 )
Inoltre, è necessario conoscere la matrice Hessiana della funzione Lagrangeana per
il problema di partenza:
2u 0
2
2
2
2
∇x,x L(x, u) = ∇x,x x1 − x2 + u(x1 + x2 − 1) =
0 2u
E’ da notare che, per ogni u > 0 l’Hessiana è definita positiva. Quindi, scelto
un punto di partenza x0 , si risolve la sequenza di problemi di programmazione
quadratica:
1
· (dk )T · ∇2x,x L(x, uk ) · dk
2
s.t. g(xk ) + (∇g(xk ))T · dk ≤ 0
minimize (∇f (xk ))T · dk +
che per l’esempio in esame ha la forma:
k
1
2u
0
k T
minimize (1, −1) · d + · (d ) ·
· dk
0 2uk
2
s.t. ((xk1 )2 + (xk2 )2 − 1) + (2x1 , 2x2 )T · dk ≤ 0
T
Al termine di ogni iterazione k si hanno a disposizione la direzione di ottimizzazione
dk ed il moltiplicatore del vincolo di disuguaglianza ūk . Si pone uk+1 = ūk ed
xk+1 = xk + τ k dk (dove τ k è un parametro che determina la lunghezza del passo), e
si itera.
Documento preparato da E. Amaldi e A. Ceselli
11
Complementi di Ricerca Operativa
Esercitazione 6
Prof. E. Amaldi
Nell’esercizio in esame (ponendo τ k = 1.0) si ottengono i seguenti risultati:
iter.
0
1
2
3
10
u
d
val.
x
(−1, 0)
1
(0, 0.25)
−1.0
(−1, 0.25)
0.5
(0.1176, 0.3456) −1.25
(−0.8824, 0.5956) 0.6176 (0.1129, 0.0554) −1.4779
(−0.7695, 0.6509) 0.7247 (0.0322, 0.0259) −1.4204
...
(−0.7076, 0.7066) 0.7071 < (10−3 , 10−3 ) −1.4142
Documento preparato da E. Amaldi e A. Ceselli
12
Esercitazione 6
Complementi di Ricerca Operativa
Prof. E. Amaldi
6.6 La funzione fmincon di MatLab permette di risolvere problemi di ottimizzazione
non lineare vincolata, utilizzando un metodo tipo programmazione quadratica sequenziale. Nel caso del problema in esame, la funzione obiettivo può essere definita
come segue (file ’cro obj.m’):
function f = cro_obj(x)
f = exp(x(1) * x(2) * x(3) *x(4) * x(5))
- 0.5 * ( (x(1))^3 + (x(2))^3 + 1)^2;
Inoltre, è necessario preparare una funzione che restituisce due vettori: un vettore
(c) contenente gli slack rispetto ai vincoli di disuguaglianza ed un vettore (ceq)
contenente l’entità della violazione dei vincoli di uguaglianza (file ’cro cons.m’, ’x’
rappresenta il vettore delle variabili, mentre ’rhs’ è il vettore dei termini noti nei
vincoli di uguaglianza):
function [c,ceq] = cro_cons(x,rhs)
c = [];
ceq = [
(x(1))^2 + (x(2))^2 + (x(3))^2 + (x(4))^2 + (x(5))^2 - rhs(1),
x(2) * x(3) - 5 * x(4) * x(5) - rhs(2),
(x(1))^3 + (x(2))^3 - rhs(3)
];
La funzione fmincon accetta diversi parametri in ingresso (utilizzare ’help fmincon’
per una descrizione completa). Ad esempio, la seguente chiamata risolve il problema
proposto nel testo, utilizzando il punto iniziale indicato; i vettori vuoti ([]) prendono il posto delle matrici dei vincoli e vettori di termini noti dei vincoli lineari,
e dei vettori che rappresentano lower e upper bound sul valore delle variabili. La
struttura ’options’ può essere creata con la seguente chiamata:
>> options = optimset(@fmincon)
e modificata come segue:
>> options = optimset(options,’Display’,’iter’)
Ad esempio, nei successivi test, è stato incrementato il livello di output della funzione
e diminuita la tolleranza sulla violazione dei vincoli.
>> [xstar,fval] = fmincon(
@(x) cro_obj(x),
[-1.8;1.7;1.9;-0.6;-0.8],
[],[],[],[],[],[],
@(x) cro_cons(x,[10;0;-1]),
options
)
Documento preparato da E. Amaldi e A. Ceselli
13
Esercitazione 6
Complementi di Ricerca Operativa
Prof. E. Amaldi
Il risultato ottenuto è il seguente:
max
Directional
First-order
Iter F-count
f(x)
constraint
Step-size
derivative
optimality Procedure
0
6
0.0580965
0.83
Infeasible start point
1
13
0.052893
0.03303
1
-0.00479
0.71
2
20
0.054456
0.0007379
1
0.00153
0.331
3
27
0.0542098
0.001425
1
-0.000167
0.304 Hessian modified
4
35
0.0540602
0.002166
0.5
-0.000138
0.301 Hessian modified
5
42
0.0539867
0.001767
1
2.41e-05
0.29 Hessian modified
6
49
0.0539499
2.36e-06
1
-3.66e-05
0.29 Hessian modified
7
56
0.0539499
9.241e-08
1
-5.8e-08
0.29 Hessian modified
Optimization terminated: Magnitude of directional derivative in search
direction less than 2*options.TolFun and maximum constraint violation
is less than options.TolCon.
xstar =
-1.7169
1.5954
1.8278
-0.7637
-0.7637
fval =
0.0539
Si nota anche come, ad esempio utilizzando un diverso punto iniziale, si possano
ottenere risultati sostanzialmente diversi:
>> [xstar, fval] = fmincon(@(x) cro_obj(x),[-1.8;1.7;1.9;0;0],[],[],[],[],[],[],
@(x) cro_cons(x,[10;0;-1]),options)
Warning: Large-scale (trust region) method does not currently solve this type of problem,
switching to medium-scale (line search).
> In fmincon at 260
max
Directional
First-order
Iter F-count
f(x)
constraint
Step-size
derivative
optimality Procedure
0
6
0.99672
3.23
Infeasible start point
1
13
-137.045
47.28
1
-699
1.28e+03
2
20
-12.4909
10.67
1
36
70.7
3
27
-0.136493
1.8
1
3.53
88.9 Hessian modified
4
34
0.980447
0.1978
1
0.221
17.9
5
41
0.999997
0.002574
1
0.000496
5.78 Hessian modified
6
48
1
2.905e-07
1
1.62e-09
0.0635 Hessian modified
Optimization terminated: first-order optimality measure less than options.TolFun
and maximum constraint violation is less than options.TolCon.
xstar =
-2.2692
2.2025
-0.0000
0
0
fval =
1.0000
Riferimenti bibliografici
[Bert03] Bertsekas D.P. (2003) Nonlinear programming, Athena scientific - Belmont
(USA)
Documento preparato da E. Amaldi e A. Ceselli
14