Universidade Estadual de Campinas~
FACULDADE DE ENGENHARIA ELETRICA E DE COMPUTACAO
DT-FEEC-UNICAMP
Metodo de Pontos Interiores em Programac~ao Linear
1 Introduc~ao
1.1 Pontos Interiores X Simplex
Ambos e cientes
Simplex: muitas iterac~oes "simples" pelas arestas
Pontos Interiores: poucas iterac~oes "caras" pelos pontos interiores.
simplex
otimo
pontos interiores
Figure 1: Simplex X Pontos Interiores
1.2 Notac~ao
A B C::: onde A = faij g - Matrizes
a b c::: onde a = faig - Vetores colunas
1.3 Programac~ao Linear - Formas Padr~oes
8
>> min
><
Primal : >> sa
>:
8
>> max
><
Dual : > sa
>:
ct x
Ax = b
x0
8
>> max
bty
bt y
><
Aty c ou Dual : >> sa Aty + z = c
>:
y livre
z 0 y livre
1.4 Denic~oes
Def. 1: Ponto interior: xo tal que xo > 0 e ponto interior
do primal.
Def. 2: Ponto factvel: xo tal que Axo = b xo 0 e um
ponto fact
vel do primal.
Def. 3: Ponto interior factvel (satisfaz 1 e 2): xo tal que
Axo = b xo > 0 e um ponto interior fact
vel do primal.
Def. 4: Gap (ou GAP): diferenca entre o valor do primal e
do dual. Ex.: GAP = ctx ; bty.
1.5 Condic~oes de Otimalidade
i) Primal factibilidade:
b ; Ax = 0 x 0
ii) Dual factibilidade:
c ; At y ; z = 0 z 0
iii) Condic~ao de Complementaridade: xizi = 0
1.6 Problemas de Mnimos Quadrados
Problema: Minx2R(x) = 12 k b ; Ax k2
Logo:
(x) = 21 (b ; Ax)t(b ; Ax) = 12 (btb ; btAx ; xtAtb + xtAtAx)
r(x) = 0 )
x = (AtA);1Atb
r2(x) = AtA > 0 )
e m
nimo global.
f(x)
x0
x1
x2
x
Figure 2: Metodo de Newton
1.7 Metodo de Newton
Caso monovari
avel:
(x) = (x0) + rt(x0)(x ; x0)
(x) = 0
) rt(x0)x = rt(x0)x0 ; (x0)
x = x0 ; r((xx00))
Caso Multivari
avel:
M
etodo de Newton:
(x) = (xk ) + r(xk )(x ; xk ) + 21 (x ; xk )tJ (xk )(x ; xk )
onde J (xk ) = r(r(xk )) (Jacobiana de (x))
r(x) = 0 ) r(xk ) + J (xk )(x ; xk ) = 0
xk+1 = xk ; J (xk )];1r(xk)
Lema 1: Uma direc~ao d para melhorar a func~ao objetivo e
fact
vel para restric~oes de igualdade do tipo Ax = b se Ad = 0
(pois A(x + d) = b ) Ad = 0).
Premissa 1: Algoritmos de pontos interiores comecam em
pontos interiores fact
veis e movem de ponto em ponto interior
fact
vel, em direca~o a soluc~ao otima.
2 Metodo Primal-Am-Escala
2.1 Base Teorica
2
66 x1
66 0
66
X (nxn) = diag(x) = 666 :
66 :
64
Miny z 12 k Xz k2
0 ::: 0
x2 ::: 0
: ::: :
: ::: :
0 0 ::: xn
sa z = c ; Aty
3
77
77
77
77
77
77
5
Miny z 12 k Xz k2 = Minz 21 k X (c ; Aty) k2
Z = k X (c ; Aty) k2
= ctX tXc ; ctX tXAty ; yt AX tXc + ytAX tXAty
X t = X , ry Z = 0: ;AX 2c + AX 2Aty = 0
z = c ; At(AX 2At);1AX 2c
Teorema (Dikin): Dados x tal que Ax = b x > 0 z =
c ; Aty y = (AX 2At);1AX 2c, ent~ao a direc~ao d = ;X 2z e uma
direc~ao de descida factvel.
Assim, d = ;X 2z minimiza k Xz k2:
xk+1 = xk + k dk com:
k
k = minik<0(; xiki ) 0 < < 1 dk = fik g
2.2 Algoritmo
1. Dados 2 (0 1) e x0 j Ax0 = b x0 > 0 k = 0:
2. Faca at
e convergir:
(a) yk = (A(X k )2At);1A(X k )2c
(b) z k = c ; Atyk
(c) dk = ;(X k )2z k
k
(d) k = minik <0f; xki g
i
(e) xk+1 = xk + k dk
(f) k k + 1
3. Fim Faca
a) Crit
erios de converg^encias:
{ Gap relativo: 1+kkbXtykkz+kcktxk k { Primal factibilidade: kb1+;Axkbkk k { Dual factibilidade: kc;A1+tykkc;k zk k { Variaca~o do valor da func~ao objetivo:
kct xk+1 ;ct xk k 1+kctxk k
b) Ponto inicial interior fact
vel x0, resolver:
8
>> min ct x + M
><
(PM ) : > sa Ax + p =
>:
(x ) b
0
onde: p = b ; Ax0. Iniciar com (x0 1).
c) C
alculo de yk :
(A(X k )2At)yk = A(X k )2c
d) am: espaco a m, escala: x~ = X ;1x:
8
>> min
><
(P ) >> sa
>:
ctx
Ax b
x0
c~ = Xc, A~ = AX .
!
x = X x~
8
>> min
><
~
) (P ) > sa
>>
:
c~tx~
A~x~ b
x~ 0
2.3 Projec~ao Am-Escala
Problema:
8>
Direc~ao d tal que:
>> max
<
>> sa
>:
8
>> max
><
>> sa
>:
ctx
Ax = b
x0
ct(x0 + d)
A(x0 + d) = b
k d k2 = 1
Lagrangeano:
L(d y ) = ct(x0 + d) + (1 ; dtd) ; yt(A(x0 + d) ; b)
Dual:
Mind y L(d y ).
Otimalidade:
i)
@L
@d
= c ; 2d ; Aty = 0
ii)
@L
@
= 1 ; dt d = 0
iii)
@L
@y
= A(x0 + d) ; b = 0
iii): Ax0 ; b + Ad = 0
)
Ad = 0
= 12 :
i): c ; d ; Aty = 0
)
)
d = c ; Aty
y = (AAt);1Ac
d = c ; At(AAt);1Ac = (I ; At(AAt);1A)c = Pc
P = (I ; At(AAt);1A) = matriz de projec~ao ortogonal
ao espaco nulo de A.
Figura 3: d = da + db, da = Atx
Adb = A(d ; da) = 0 (espaco nulo de A).
Assim:
Ad = Ada = AAtx
x = (AAt);1Ad
)
da = At(AAt);1Ad
db = d ; At(AAt);1Ad = (I ; At(AAt);1A)d = P (A)d
da
db
N(A)
d
t
R(A )
Figure 3: Matriz de proje
c~ao ortogonal de A
Problema (P~ ):
d~ = ;P~ c~ = ;(I ; A~t(A~A~t);1A~)~c
A~ = AX c~ = Xc, iterac~ao k:
d~k = ;(I ; X k At(AX k X k At);1AX k )X kc
= ;X k c + X k At(A(X k )2At);1A(X k )2c
Logo:
8 k
>< d
>: xk+1
= X k d~k = ;(I ; (X k )2At(A(X k )2At);1A):(X k )2c
= xk + k dk
Primal-a m-escala:
8
>> dk
><
k
>> z
: yk
)
8 k
>< d
>: xk+1
= ;(X k )2z k
= c ; Atyk
= (A(X k )2At);1A(X k )2c
= ;(I ; (X k )2At(A(X k )2At);1A):(X k )2c
= xk + k dk
2.4 Exemplo (Frannie's Firewood Problem)
Frannie vende 3 "cordas" de lenha todo nal do ano. Pode vender
meia "corda" a U $90 ou uma "corda" a U $150. Como maximizar
o lucro?
8
>> max 90x + 150x
1
2
><
1
Modelo : > sa 2 x1 + x2 3
>:
x1 x2 0
x3 = vari
avel de folga, (1 21 2) = ponto inicial interior fact
vel
Resultado: tabela 1, gura 4.
x1
x2
x3
ctx
1.00
1.73
2.82
3.92
5.66
5.88
5.98
5.99
6.00
6.00
0.50
0.80
1.13
1.00
0.14
0.05
0.01
0.005
0.001
0.000
2.00
1.34
0.46
0.04
0.03
0.01
0.003
0.000
0.000
0.000
165.00
275.51
423.40
502.84
530.45
537.25
539.22
539.79
539.95
539.99
c1 = 90 c2 = 150 c3 = 0 b = 3
550
500
450
c’x
400
350
300
250
200
150
1.5
6
1
5
4
0.5
3
2
x2
0
1
x1
Figure 4: Evolu
c~ao dos pontos interiores do problema de Frannie
3 Metodo Dual-Am-Escala
3.1 Introduc~ao
Problema:
8>
< min 12 k Zx k2
>: sa Ax = b
Z = diag(zi)
Lagrangeano:
L(x w) = 21 k Zx k2 +wt(b ; Ax) = 12 (xtZZx) + wt(b ; Ax)
Condic~oes de otimalidade:
Ou seja:
8> @L
< @x
>: @L
@w
=0
=0
!
!
Z 2 x ; At w = 0
b ; Ax = 0
x = Z ;2Atw (AZ ;2At)w = b ) x = Z ;2At(AZ ;2At);1b
Hessiana:
H
2
= 64
3
2
Z ;A 75 = ;AAt
;At 0
z = c ; Aty ) dz = ;Atdy
<0
Direca~o dz = ;Z 2x:
dz = ;Z 2Z ;2At(AZ ;2At);1b]
= ;At(AZ ;2At);1b = ;Atdy
Assim:
dy = (AZ ;2At);1b
Teorema: Dados (y z) tais que z = c ; Aty, z 0, Ax = b,
x = Z ;2At(AZ ;2At);1b, a direc~ao dada por:
(dy dz ) = (AZ ;2At);1b ;Z 2x)
e dual fact
vel e e de subida.
3.2 Algoritmo
1. Dados (y0 z 0) tal que Aty0 + z 0 = c z 0 > 0 e 0,
2. Faca at
e converg^encia:
(a) dyk = (A(Z k );2At)b
(b) dz k = ;Atdyk
(c) xk = ;(Z k );2dz k
2 (0
1) k =
k
(d) k = minzik <0f; zzik g
i
(e) yk+1 = yk + k dyk
(f) z k+1 = z k + k dz k (ou z k+1 = c ; Atyk+1)
(g) k k + 1
3. Fim Faca
a) Crit
erio de converg^encia originalmente utilizado:
kbt y k ;bt y k+1k .
max(1 kbt yk k)
b) Ponto inicial dual fact
vel (y0 z 0), resolver o problema:
8
>> max
><
>> sa
>:
bty ; M
Aty + z ; e = c
z0
aplicando o m
etodo dual-a m-escala at
e < 0, com valor
inicial y0 qualquer (livre), 0 = ;2minj (cj ; Atj y0).
c) z k+1 = c ; Atyk+1 = c ; At(yk + k dyk ) =
c ; Atyk ; k Atdyk = zk + k dzk
d) O c
alculo de xk e dispens
avel, a n~ao ser que se utilize no
crit
erio de converg^encia.
e) Como y e livre, n~ao e feito teste de barreira.
f) Grande custo computacional no c
alculo de A(Z k );2At.
4 Metodo Primal-Dual-Am-Escala
4.1 Introduc~ao
Condic~oes de otimalidade:
F (x y
2
66 Fp
z) = 6664 Fd
Fa
3 2
3
Ax
;
b
77 66
7
77 = 66 At y + z ; c 777 = 0
75 64
75
Aproximac~ao linear fornece:
XZe
(x y z ) = (x0 y0 z 0) ; J ;1(x0 y0 z 0)F (x0 y0 z 0)
pois:
F (x y z) = F (x0 y0 z 0) + J (x0 y0 z 0)(x y z );
(x0 y0 z 0))] = 0
e:
2
66 b ; Ax0
;F (x0 y 0 z 0 ) = 666 c ; At y 0 ; z 0
4
;X 0Z 0e
2
t
66 rFp
J (x0 y0 z0) = 6664 rFdt
rFat
3 2
77 66
77 = 66
75 64
3 2
77 66 rp
77 = 66 r
75 64 d
ra
3
77
77 = r
75
3
A 0 0 77
0 At I 7775
Z0 0 X0
Assim:
(x y
2
66 A
0
0
0
z) = (x y z ) + 6664 0
Z0
3;1 2
0 0 77 66 rp
At I 7775 6664 rd
0 X0
ra
3
77
77 = (x0 y 0 z 0) + d
75
onde:
2
3 2
66 dx 77 66 A
d = 6664 dy 7775 = 6664 0
dz
Z0
3;1 2
0 0 77 66 rp
At I 7775 6664 rd
0 X0
ra
Pode-se resolver o sistema:
2
66
66
64
32
3 2
A 0 0 77 66 dx 77 66 rp
0 At I 7775 6664 dy 7775 = 6664 rd
Z 0 0 X 0 dz
ra
3
77
77
75
Considerando:
8
>> Adx
><
t
>> A dy + dz
>: Z 0dx + X 0dz
dz = rd ; Atdy
)
= rp
= rd
= ra
:
Z 0dx + X 0dz = Z 0dx + X 0(rd ; Atdy) = ra
3
77
77
75
Z 0dx ; X 0Atdy = ra ; X 0rd
Ou seja:
;(X 0);1Z 0 dx + Atdy
= ;(X 0);1ra + rd
que fornece:
2
64
32
3
2
3
A 0 75 64 dx 75 = 64 rp
75
t
0
;
1
;D A
dy
rd ; (X ) ra
onde D = X ;1Z . Note que:
8
>><> Adx
t
>> ;Ddx + A dy
>:
e:
dx
= rp
= rd ; X ;1ra
= D;1Atdy ; rd + X ;1ra]
D;1Atdy = dx + D;1(rd ; X ;1ra)
(AD;1At)dy = Adx + AD;1rd ; AD;1X ;1ra
)
(AD;1At)dy = rp + AD;1rd ; AZ ;1ra
4.2 Algoritmo
1. Dados (x0 y0 z 0) tal que (x0 z 0) > 0 e 2 (0 1) k = 0,
2. Faca at
e converg^encia:
(a) rpk = b ; Axk
(b) rdk = c ; Atyk ; z k
(c) rak = ;X k Z k e
(d) Dk = (X k );1Z k ]
(e) dyk = A(Dk );1At];1rpk + A(Dk );1rdk ; A(Z k );1rak ]
(f) dxk = (Dk );1Atdyk ; rdk + (X k );1rak ]
(g) dz k = (X k );1rak ; Z k dxk ]
xki
k
(h) p = minxki <0f; xk g
i
k
minzik <0f; zziik g
(i) kd =
(j) pk = min(
kp 1)
(k) dk = min(
kd 1)
(l) xk+1 = xk + pk dxk
(m) yk+1 = yk + dk dyk
(n) z k+1 = z k + dk dz k
(o) k k + 1
3. Fim Faca
A converg^encia pode ser testada sobre o valor de k F k2.
O ponto inicial (x0 y0 z 0) n~ao precisa ser fact
vel. Recomendase:
{ Para o primal:
Fazendo x = Aty~, Ax = b
!
AAty~ = b
) y~ = (AAt );1b
Aty~ = At(AAt);1b
x = At(AAt);1b
1 = max(;minixi 2 2kkbAk1k1 ), onde 2 = 100
x0i = max(xi 1)
(k b k1 = Pmi=1 j bi j k A k1 = maxj Pmi=1 j aij j)
{ Para dual:
Aty + z = c y livre e z > 0
8
>
>
>
< zi + 3
0
z = >> ;zi
>
: 3
)
y0 = 0, e:
se zi 0
se zi ;3
se ;3 zi 0
onde 3 = 1 + k c k1. Estes pontos procuram ser bem
posicionados, longe da fronteira (xizi n~ao muito pequenos).
O tamanho do passo para y e o mesmo de para z (dk ) para
garantir que ocorra c ; Atyk ; z k = 0 na converg^encia:
c ; Atdyk+1 ; zk+1 = c ; At(yk + dk dyk ) ; (zk + dk dzk )
= (c ; Atyk ; z k ) ; dk (Atdyk ; dz k )
Como n~ao precisa de ponto inicial fact
vel, este m
etodo e melhor que o primal ou dual (n~ao precisa de fase I).
5 Metodo Primal-Dual Classico
5.1 Introduc~ao
Primal-dual a m-escala permite que x e z aproximem rapidamente
das fronteiras ) ine ciente.
Primal-Dual Cl
assico acrescenta uma perturbac~ao na condica~o
de complementariedade:
xi zi = i
Novas condico~es de otimalidade:
8
>> b ; Ax
><
t
t
>> c ; A y ; z
>: e ; XZe
= 0
= 0
=0
tal que limk ! 1k = 0
Estimac~ao de k :
k = Tr(X k )tZ k ] onde: TrX ] = traco de X
Na maioria das implementac~oes, adota-se:
k = k ( nk ) k 2 (0 1)
Nota: k = 0
)
primal-dual a m-escala, k = 1
)
direca~o de centragem pois:
e ; XZe = 0
k e ; XZe = 0
n
X tZ e
n
A cada iterac~ao, tem-se o sistema J (xk yk z k )dk = rk , ou seja:
2
66
66
64
32 k
A 0 0 77 66 dx
0 At I 7775 6664 dyk
Z k 0 X k dzk
) X tZe =
3 2 k
77 66 rp
77 = 66 rk
75 64 d
rck
3 2
3
k
77 66 b ; Ax
77
77 = 66 c ; At ; z k
77
75 64
75
k
k
t
k
e ; (X ) Z e
Tem-se assim, duas diferencas em relac~ao ao m
etodo primal-dual
a m-escala:
a) troca de rak por rck
b) c
alculo de k
5.2 Algoritmo
1. Dados (x0 y0 z 0) tal que (x0 z 0) > 0, 2 (0 1) 2 (0 1) e
k = 0,
2. Faca at
e converg^encia:
(a) k = TrX k Z k ]
k
(b) k = ( n )
(c) rpk = b ; Axk
(d) rdk = c ; Atyk ; z k
(e) rck = k e ; X k Z k e
(f) Dk = (X k );1Z k ]
(g) dyk = A(Dk );1At];1rpk + A(Dk );1rdk ; A(Z k );1rck ]
(h) dxk = (Dk );1Atdyk ; rdk + (X k );1rck ]
(i) dz k = (X k );1rck ; Z k dxk ]
xki
k
(j) p = minxki <0f; xk g
i
k
minzik <0f; zziik g
(k) kd =
(l) pk = min(
kp 1)
(m) dk = min(
kd 1)
(n) xk+1 = xk + pk dxk
(o) yk+1 = yk + dk dyk
(p) z k+1 = z k + dk dz k
(q) k k + 1
3. Fim Faca
O crit
erio de converg^encia e a inicializac~ao deste m
etodo podem ser os mesmos do m
etodo primal-dual a m-escala. E
recomend
avel inicializar com 0 alto.
Dependendo dos valores de e , obtemos algoritmos de diversas naturezas (complexidade polinomial, converg^encia superlinear, etc.).
Os valores t
picos de est~ao entre (0:995 0:99995).
k2
Quando k < 1, recomenda-se utilizar k = (n) para procurar acelerar a converg^encia.
6 Metodo Preditor-Corretor
6.1 Introduc~ao
Baseado em 3 componentes:
direc~ao a m-escala d~ (direc~ao de Newton, preditor).
direc~ao de centragem, de nido pelo do primal-dual cl
assico.
direc~ao de correc~ao d^, que tenta compensar a aproximac~ao
linear de Newton.
Id
eia: calcular a direc~ao a m-escala e estudar o progresso ao
longo desta direc~ao, atuando na perturbac~ao (centragem) e na
correca~o n~ao-linear.
No ponto (x y z ):
8
>>
< t Adx~
(1) >> A dy~ + dz~
>: Zdx~ + Xdz~
= rp
= rd
= ra = ;XZe
Obt
em-se, ent~ao, o ponto (~x y~ z~), onde:
2 3 2
3
66 x~ 77 66 x + dx~ 77
66 y~ 77 = 66 y + dy~ 77
75
64 75 64
z~
z + dz~
A seguir, determinar a direc~ao (dx^ dy^ dz^):
8
>>
Adx^
><
t
(2) > A dy^ + dz^
>: Zdx^ + Xdz^
= 0
= 0
= e ; (Dx~Dz~)e = rc
onde Dx~ = diag(dx~) e Dz~ = diag(dz~). Finalmente, a direc~ao
nal (dx dy dz ) e determinada somando (1) e (2):
8>
>>
= rp
A(dx~ + dx^)
< t
>> A (dy~ + dy^) + (dz~ + dz^) = rd
>: Z (dx~ + dx^) + X (dz~ + dz^) = ra + rc = rs
onde:
8
>< ra
>: r
c
= ;XZe
= e ; (Dx~Dz~)e
6.2 Algoritmo
1. Dados (x0 y0 z 0) tal que (x0 z 0) > 0, 2 (0 1) e k = 0,
2. Faca at
e converg^encia:
(a) rpk = b ; Axk
(b) rdk = c ; Atyk ; z k
(c) rak = ;X k Z k e
(d) Dk = (X k );1Z k ]
(e) dy~k = A(Dk );1At];1rpk + A(Dk );1rdk ; A(Z k );1rak ]
(f) dx~k = (Dk );1Atdy~k ; rdk + (X k );1rak ]
(g) dz~k = (X k );1rak ; Z k dx~k ]
k
(h) ~kp = minx~ki <0f; xx~ik g
i
k
minz~ik <0f; zz~iik g
(i) ~kd =
(j) ~ pk = min( ~kp 1)
(k) ~ dk = min( ~kd 1)
(l) ~ k = (xk + ~ pk dx~k )t(z k + ~ dk dz~k ), k = TrX k Z k ]
8 k
>> ~ 3
k>1
)
se
<
(
k
(m) k = >> k
: ( pn ) se k 1
k
(n) k = k ( n )
(o) rsk = rak + k e ; (Dx~k )(Dz~k )e
(p) dyk = A(Dk );1At];1rpk + A(Dk );1rdk ; A(Z k );1rsk ]
(q) dxk = (Dk );1Atdyk ; rdk + (X k );1rsk ]
(r) dz k = (X k );1rsk ; Z k dxk ]
xki
k
(s) p = minxki <0f; xk g
i
k
minzik <0f; zziik g
(t) kd =
(u) pk = min(
kp 1)
(v) dk = min(
kd 1)
(w) xk+1 = xk + pk dxk
(x) yk+1 = yk + dk dyk
(y) z k+1 = z k + dk dz k
(z) k k + 1
3. Fim Faca
Nota:
O crit
erio de converg^encia e a inicializac~ao deste m
etodo podem ser os mesmos do m
etodo primal-dual a m-escala.
Dois sistemas lineares precisam ser resolvidos, utilizando a
mesma relac~ao: A(Dk );1At = Lk (Lk )t.
Espera-se que o esf^orco para resolver dois sistemas lineares seja
recompensado pela reduc~ao no n
umero de iterac~oes.
Este e o m
etodo com melhores resultados te
oricos e pr
aticos
(tem converg^encia quadr
atica).
7 Metodo de Barreira Logartimica
7.1 Introduc~ao
Seja o problema:
8
>> max ct x
><
(P 1) : >> sa Ax = b
>:
x0
Substituindo a restric~ao x 0 na forma:
8
>> max ct x + tf (x)
>>
>> sa Ax = b
2
3
>>
>>
66 ln(x1) 77
>>
66 ln(x ) 77
<
2 7
66
(P~ 1) : >>
66 : 777
>> onde : f (x) = ln(x) = 66
7
>>
66 : 777
>>
66
7
>>
66 : 777
>>
4
:
ln(xn) 5
pode ser assumido um escalar ( 2 R).
Exemplo:
8
>> min
>>
<
>>
>>
:
z = 5x1 + 3x2
sa x1 + x2 1 (eq:1)
0 x1 2 (eq:2)
0 x2 2 (eq:3)
x2
eq. 2
00
11
000
0
1
0000
1111
00000
11111
00000000
11111111
2 111
00
11
000
111
0000
1111
0
1
00000
11111
00
11
0
1
000
111
0000
1111
00000
11111
eq. 1 111
00
11
0eq. 3
1
0000
1111
00000
11111
0000
1
00000
11111
000
111
0000
1111
00000
11111
0000
1
00000
11111
0000
1111
00000
11111
00000
11111
1 111
000
111
0000
1111
00000
11111
00000
11111
000
111
0000
1111
00000
11111
00000
11111
000
111
00000
11111
000
111
00000
11111
000
111
000
111
00000
11111
000
111
000
111
1
x1
2
Figure 5: Regi~ao factvel do problema
Regi~ao fact
vel na gura 5. Com as vari
aveis de folga:
8>
>> min z = 5x1 + 3x2
>>
>< sa x1 + x2 ; x3 = 1
x1 + x4 = 2
>
>>
x2 + x5 = 2
>>
>:
xi 0 i = 1 2 3 4 5
Com a adic~ao da func~ao barreira:
8
>
min z~
>
>
>
>
>
<
> sa
>>
>>
>:
= 5x1 + 3x2 ; (ln(x1) + ln(x2) + ln(x3)+
ln(x4) + ln(x5))
x1 + x2 ; x3 = 1
x1 + x4 = 2
x2 + x5 = 2
Vamos analisar duas situaco~es:
a) x1 = x2 = 54 (aproximadamente no meio da regi~ao fact
vel)
z = 5x1 + 3x2 = 10
z~ = 7:237
b) x1 = 1:999 x2 = 0:001 (quase na fronteira)
z = 9:998
z~ = 134:3
(P 1) pode ser colocado como sendo:
8>
<
(P 2) : >: max
sa
f (x) = ctx + lnx
Ax = b
e o problema de busca de melhor direc~ao fact
vel (em torno do
ponto xk ) pode ser colocado como:
8
>< max f (xk dxk ) rf t(xk )dxk + 21 (dxk )tJ (xk )dxk
=
(P 3) : >: sa Adxk = 0
Note que:
8>
>> rf (xk ) = c + x1k
>><
J (xk ) = f8aij g
>< ; se i = j
>>
(xki )2
>> aij
=
>: 0
>:
se c:c:
Aplicando o m
etodo primal a m-escala ao P (3), obtemos (do
lagrangeano associado):
dxk = 1 X k P k(ck + e)
onde X k = diag(xk ) ck = cX k e P k e a matriz de projec~ao
P k = I ; (Ak )tAk (Ak )t];1Ak com Ak = AX k .
Pr
oximo ponto interior:
xk+1 = xk + k dxk , com dxk = X k P k(ck + e)
Melhor tamanho do passo no m
etodo de Newton = 1 , ):
k
k = minf 1 maxg, onde: max = minxk <0f; xxik g
i
in!uencia fortemente no fator de converg^encia.
Valor t
pico: 2 :9
:999]
7.2 Algoritmo
1. Dados 2 (0 1), x0 > 0, 0 grande, k = 0:
2. Faca at
e convergir:
(a) X k = diagfxk g
(b) ck = X k c
(c) Ak = AX k
(d) P k = I ; (Ak )tAk (Ak )t];1Ak
(e) dxk = X k P k (ck + k e)
xki
k
(f) max = minxk <0 f; k g
i
xi
i
(g) k = minf 1k kmaxg
(h) xk+1 = xk + k dxk
(i) k+1 = f (k ctxk )
(j) k k + 1
3. Fim Faca
Nota:
a) Crit
erio de converg^encias: pode ser feito sobre a variac~ao do
valor de ctxk ( ctxk ) ou sobre o valor de dxk ( dxk ).
b) O ponto inicial interior x0 pode ser encontrado utilizando o
m
etodo de M grande, como no caso do primal a m-escala.
c) O valor de k pode ser tal que:
8
>< 0
k
= >: k
ctxk k
se
se
k
k
ctxk k grande
ctxk k pequeno
ou considerando:
k+1 = k 0 < < 1
O valor de 2 R < 1 de ne a velocidade de converg^encia.
8 Comentarios sobre Sistemas Lineares
Nos m
etodos de pontos interiores, s~ao precisos determinar matrizes
do tipo: B = AD;1At, onde A e uma matriz dada e D > 0 e uma
matriz diagonal que varia a cada iterac~ao e e tal que:
Primal A m-Escala: D = X ;2 (yk = A(X k )2At];1A(X k )2c)
Dual A m-Escala: D = Z 2 (dyk = A(Z k );2At]b)
Primal-Dual A m-Escala:
D = X ;1Z
(dyk = A(Dk );1At];1rpk + A(Dk );1rdk ; A(Z k );1rck ])
E
uma etapa computacionalmente "cara". A seguir, vamos analisar alguns aspectos que poder~ao melhorar o comportamento computacional.
Como D > 0 e diagonal, Dt > 0 tamb
em e diagonal e podemos
escrever:
B = AD;1At = A~A~t > 0 com A~ = AD; 12
Assim, o elemento ij de B e dado por (ver gura 6):
ij = (~ai)t(~aj ) = Pnk=1 a~ik a~jk = Pni=1 aik ajk kk;1
onde fij g = D. Para cada ij , temos (2 produtos + uma
soma)n] = 3n operac~oes. Como aik ajk e constante para todas as
i
j
~t
a
i
a
j
Figure 6: (~a )(~a )
t
i
j
iterac~oes, se armazenarmos este produto, o c
alculo de ij vai exigir
(1 produto + 1 soma)n] = 2n operac~oes.
Existem m(m2+1) elementos diferentes em B (sim
etrica, portanto,
m + (m ; 1) + ::: + 1 = m2+1 m). Assim, o n
umero de operac~oes
s~ao:
Sem armazenamento: 3n m(m2+1) operac~oes
Com armazenamento: nm(m + 1) operac~oes
A matriz B e sim
etrica e de nida positiva. Logo, a resoluc~ao do
sistema:
Bd = b
ca facilitada, se decompormos na forma:
LUd = b
8
>< L =
>: U =
matriz triangular inferior unitario
onde:
matriz triangular superior unitario
Aqui, unit
ario signi ca matriz com 1's na diagonal principal. A
resoluc~ao e processada, fazendo:
=
L
w
U
b
d
w
=
Figure 7: Decomposi
ca~o LU - Solu
c~ao por substitui
ca~o
L(Ud) = b
8>
< Lw = b
>: Ud = w
A decomposica~o LU n~ao e u
nica mas a decomposic~ao B = L U
e u
nica (onde e diagonal e L e U s~ao matrizes unit
arios. E
poss
vel mostrar que, para matrizes sim
etricas, U = Lt, ou seja,
B = L1 L1 t. Como B > 0 ) > 0, podemos escrever que
= 2 2 , o que permite escrever:
B=L
1
2
L = L~ L~ t
1
2 t
que e conhecido como decomposic~ao de Cholesky. Isso permite
escrevermos:
8
>< L
~w = b
t
~
~
Bd = b ) LL d = b ) >: L~ td = w
FEEC, 09 de marco de 2000.
Akebo Yamakami - Orientador
DT-FEEC-UNICAMP