Upload
vantuyen
View
220
Download
0
Embed Size (px)
Citation preview
Planejamento do Tratamento por Radioterapia
Atraves de Metodos de Pontos Interiores 1
Cecılia Bollini Barboza Cid
Orientador: Prof. Dr. Aurelio Ribeiro Leite de Oliveira
Instituto de Ciencias Matematicas e de Computacao
ICMC-USP - Sao Carlos
1Trabalho realizado com o apoio da CNPq
2
Agradecimentos
Quero manifestar minha gratidao ao meu orientador Professor Dr. Aurelio Ribeiro
Leite de Oliveira pela amizade, crıticas, sugestoes e apoio, os quais proporcionaram a
elaboracao deste trabalho.
Meus agradecimentos a minha famılia, especialmente aos meus avos, Plınio e Izelda,
por todo seu amor, sabedoria, incentivo e oracoes.
Um obrigado muito especial ao Alexandre que me apoiou em todo momento, com
muito amor, compreensao e incentivo.
Tambem gostaria de agradecer a todos meus companheiros de mestrado que con-
tribuıram para realizacao deste trabalho e que se revelaram grandes amigos.
Quero manifestar meu reconhecimento e gratidao a Deus, que sempre me mostrou
o caminho a ser seguido, nunca deixou que eu me esquecesse do quanto a vida e realmente
bela, e que e verdadeiramente o responsavel pela realizacao deste trabalho.
Finalmente, gostaria de agradecer ao apoio financeiro do CNPq e a todos que de
diversas formas contribuıram para a realizacao deste projeto de mestrado.
Resumo
O objetivo deste trabalho consiste no desenvolvimento, estudo e implementacao de metodos
de pontos interiores especıficos para o problema de planejamento do tratamento de cancer
por radioterapia. Este e um problema de grande porte que contem uma estrutura matri-
cial particular. A exploracao desta estrutura de forma eficiente obtem bom desempenho
computacional, atraves da reducao da dimensao dos sistemas lineares que devem ser re-
solvidos a cada iteracao, agilizando a definicao de um tratamento adequado, uma vez que
tipicamente varias simulacoes sao realizadas antes da definicao de um plano definitivo.
Resultados numericos em Matlab ilustram a eficiencia desta abordagem em problemas
reais e mostram a superioridade do metodo preditor corretor em comparacao ao metodo
primal-dual.
i
ii
Abstract
In this work, a specialized interior point method is developed for planning cancer treat-
ment by radiotherapy. This is a large-scale problem with a specific matrix structure. That
structure is explored in an efficient way reducing the dimension of the linear system which
must be solved at each iteration speeding up the treatment design since usually several
versions must be solved to obtain a satisfactory plan. Moreover, the system obtained is
sparse, symmetric and positive definite. Numerical results in Matlab illustrate the effi-
ciency of this approach in real problems and show the superiority of the preditor-corrector
method in comparison to the primal-dual method.
iii
iv
Conteudo
1 Introducao 1
1.1 Abordagem por Metodos de Pontos Interiores . . . . . . . . . . . . . . . . 1
2 Metodos de Pontos Interiores para Programacao Linear 3
2.1 Metodo de Newton para uma Variavel . . . . . . . . . . . . . . . . . . . . 4
2.2 Metodo de Newton para Varias Variaveis . . . . . . . . . . . . . . . . . . . 5
2.3 Metodo Primal-Dual Afim Escala . . . . . . . . . . . . . . . . . . . . . . . 5
2.3.1 Metodo Primal-Dual Classico . . . . . . . . . . . . . . . . . . . . . 7
2.4 Metodo Preditor-Corretor . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.5 Calculo das Direcoes nos Metodos de Pontos Interiores . . . . . . . . . . . 10
3 Um Modelo para Formulacao do Tratamento por Radioterapia 13
3.1 Formulacao do Modelo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.2 Problema de Otimizacao Linear . . . . . . . . . . . . . . . . . . . . . . . . 18
3.3 Propriedades das Restricoes Elasticas . . . . . . . . . . . . . . . . . . . . . 20
3.3.1 Interpretacao das Solucoes . . . . . . . . . . . . . . . . . . . . . . . 20
4 Metodos de Pontos Interiores Aplicados ao Modelo de Tratamento por
Radioterapia 23
4.1 Analise Media . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4.1.1 O Metodo Primal-Dual para Analise Media . . . . . . . . . . . . . . 28
4.1.2 Eliminacao de Variaveis . . . . . . . . . . . . . . . . . . . . . . . . 30
4.1.3 Resumo do Metodo de Pontos Interiores . . . . . . . . . . . . . . . 35
4.1.4 Metodo Preditor-Corretor para Analise Media . . . . . . . . . . . . 36
4.2 Analise Absoluta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
v
4.2.1 O Metodo Primal-Dual para Analise Absoluta . . . . . . . . . . . . 43
4.2.2 Eliminacao de Variaveis . . . . . . . . . . . . . . . . . . . . . . . . 45
4.2.3 Resumo do Metodo de Pontos Interiores . . . . . . . . . . . . . . . 52
4.2.4 Metodo Preditor-Corretor para Analise Absoluta . . . . . . . . . . 53
5 Resultados Computacionais 55
6 Conclusoes 59
6.1 Perspectivas Futuras . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
Bibliografia 61
vi
Capıtulo 1
Introducao
Modelos de programacao linear tem sido usados extensivamente para encontrar bons
planos de tratamento por radioterapia. A primeira formulacao proposta como um modelo
de programacao linear ocorreu em 1968 [3]. O objetivo do tratamento e eliminar as celulas
do tumor atraves de radiacao ao mesmo tempo em que procura evitar a destruicao de
celulas vizinhas saudaveis tambem afetadas pela radiacao. Existem atualmente maquinas
sofisticadas inspiradas na tomografia computadorizada que emitem radiacao ao longo de
todo o corpo do paciente [23], por um lado ampliando as possibilidades de minimizacao
destes efeitos colaterais e por outro aumentando a complexidade de obtencao de planos
de tratamento.
Do ponto de vista matematico, o desafio consiste em emitir uma alta dosagem
de radiacao no tumor, suficiente para sua eliminacao, e simultaneamente, minimizar a
radiacao nas regioes vizinhas compostas de tecido saudavel, reduzindo ao maximo as
complicacoes nestas regioes que sao muitas vezes crıticas.
1.1 Abordagem por Metodos de Pontos Interiores
Uma vantagem deste problema, do ponto de vista computacional, e que ele pode ser mo-
delado como um problema de otimizacao linear [3, 8, 9, 22, 24, 23]. Por outro lado, o tipo
de solucao obtida pelo metodo simplex para estes modelos nao corresponde ao tratamento
mais adequado do ponto de vista medico porque as variaveis nao basicas estarao nos seus
respectivos limites [8]. Isto significa, nestes modelos, que parte do tumor recebe a dosagem
1
mınima estabelecida para sua eliminacao enquanto parte do tecido saudavel pode receber
a dosagem maxima permitida.
Solucoes deste tipo sao naturalmente evitadas pelos metodos de pontos interiores.
Na presenca de multiplas solucoes esses metodos convergem para aquelas, em certa me-
dida, distanciadas de vertices. Estas seriam as solucoes mais adequadas na formulacao
de um tratamento, pois as variaveis nao tendem a atingir seus limites na presenca de
multiplas solucoes. Uma outra vantagem da abordagem por pontos interiores e a possi-
bilidade de trabalhar com pontos infactıveis durante o processo de otimizacao, facilitando
a identificacao de modelos mal formulados.
Uma motivacao adicional para a abordagem por metodos de pontos interiores,
esta na estrutura matricial bem definida deste problema. Experiencias anteriores com o
desenvolvimento de metodos de pontos interiores especıficos para uma classe de proble-
mas e a exploracao da estrutura matricial dos sistemas lineares resultantes mostra que
esta abordagem e muito bem sucedida em termos da obtencao de melhor desempenho
computacional [4, 15, 16, 18, 17, 19, 20, 21].
Este trabalho esta organizado da seguinte forma: no Capitulo 2 sao discutidos
os metodos de pontos interiores para programacao linear. No Capitulo 3 apresentamos
modelos para a formulacao do tratamento por radioterapia e discutimos com mais detalhes
o modelo adotado. No Capıtulo 4 os metodos de pontos interiores sao aplicados ao
modelo de tratamento por radioterapia. No Capitulo 5 sao apresentados os resultados
computacionais obtidos. Finalmente, no Capitulo 6 sao apresentadas as conclusoes e
perspectivas futuras deste trabalho.
2
Capıtulo 2
Metodos de Pontos Interiores para
Programacao Linear
Neste capıtulo vamos adotar a seguinte forma, denominada padrao, para o problema de
programacao linear:
min cT x
s.a Ax = b
x ≥ 0,
onde, A ∈ <m×n com posto m e os vetores coluna x, c e b tem a dimensao apropriada.
O dual desse problema e dado por:
max bT y
s.a AT y + z = c
z ≥ 0.
Um ponto e dito interior quando todas as variaveis estao estritamente dentro de
seus limites [27] no problema dual por exemplo z0 > 0 e um ponto interior. Um ponto
interior e factıvel quando todas as restricoes sao satisfeitas, ou seja, o ponto x0, onde,
Ax0 = b com x0 > 0 e um ponto interior factıvel. As condicoes de otimalidade para
problemas primal e dual sao dadas pela [27]:
3
factibilidade primal
b − Ax = 0
x ≥ 0,
factibilidade dual
c − AT y − z = 0
z ≥ 0,
e condicoes de complementaridade [26]
xizi = 0, i = 1 · · ·n
ou seja, XZe = 0,
onde, X e a matriz diagonal formada pelos elementos do vetor x, Z e a matriz diagonal
formada pelos elementos do vetor z e e e o vetor unitario, eT = (1, 1, 1, · · · , 1). Os metodos
de pontos interiores primais-duais podem ser desenvolvidos atraves aplicacao do metodo
de Newton as condicoes de otimalidade [27] partindo de um ponto interior desconsiderando
as desigualdades. A nao negatividade das variaveis e controlada pelo tamanho do passo.
2.1 Metodo de Newton para uma Variavel
Antes de desenvolvermos estes metodos, apresentaremos uma breve revisao do metodo de
Newton. O metodo de Newton para uma variavel busca zeros de uma funcao resolvendo
equacoes da forma, φ(x) = 0. Este metodo pode ser deduzido aplicando a serie de Taylor
a φ(x) em torno de x0, obtendo:
φ(x) = φ(x0) + φ′(x0)(x − x0) + · · ·
Ignorando os termos de ordem superior temos, para φ(x) = 0 com φ′(x) 6= 0 :
− φ(x0)φ′(x0)
= x − x0 ⇒ x = − φ(x0)φ′(x0)
+ x0
obtendo assim o processo iterativo:
xk+1 = xk − φ(xk)φ′(xk)
, onde dk = − φ(xk)φ′(xk)
e a direcao de Newton.
Este processo pode ser repetido ate que uma tolerancia estabelecida seja satisfeita.
4
2.2 Metodo de Newton para Varias Variaveis
O metodo de Newton para varias variaveis tem como objetivo encontrar os zeros de
sistemas de equacoes da seguinte forma:
seja, f(x)=0, um conjunto de funcoes nao lineares φi(x) onde φi(x) = 0, i = 1, · · · , n.
O metodo pode ser deduzido aplicando a serie de Taylor para cada φi(x) em torno de x0,
obtendo:
φi(x) = φi(x0) + [∇φi(x
0)]T (x − x0) + · · · , para todo, i = 1 · · ·n, onde
∇φi(x0) =
∂φi(x0)
∂x1
...
∂φi(x0)
∂xn
Ignorando os termos de ordem superior temos:
[∇φi(x0)]T (x − x0) = −φi(x
0), i = 1, · · · , n.
Onde definimos agora o Jacobiano no ponto x0:
J(x0) =
∇φ1(x0)
...
∇φn(x0)
T
e f(x0) =
φ1(x0)
...
φn(x0)
Podemos resumir assim o metodo de Newton:
Dado um ponto inicial x0
Para k = 0, 1, · · ·
Calculo da direcao dk = −J(xk)−1f(xk)
Calculo do novo ponto xk+1 = xk + dk
ate convergir
onde um criterio de convergencia poderia ser ‖f(xk+1)‖/‖xk+1‖ menor que uma dada
tolerancia.
2.3 Metodo Primal-Dual Afim Escala
Este metodo busca a solucao dos problemas de programacao linear primal e dual aplicando
variantes do metodo de Newton as condicoes de otimalidade, e modificando o tamanho
do passo para manter os pontos interiores [14]. Escrevemos as condicoes de otimalidade
5
da seguinte forma,
f(x, y, z) =
Ax − b
AT y + z − c
XZe
= 0.
Aplicando o metodo de Newton a este sistema de equacoes nao lineares obtemos:
(x, y, z) = (x0, y0, z0) − (J(x0, y0, z0))−1f(x0, y0, z0) = (x0, y0, z0) − d0.
onde,
−f(x0, y0, z0) = −
Ax0 − b
AT y0 + z0 − c
X0Z0e
=
r0p
r0d
r0a
≡ r0.
A direcao do metodo de Newton aplicado as equacoes f(x,y,z) torna-se:
J(x0, y0, z0)d0 =
A 0 0
0 AT I
Z0 0 X0
d0x
d0y
d0z
=
r0p
r0d
r0a
logo,
d0 =
d0x
d0y
d0z
=
A 0 0
0 AT I
Z0 0 X0
−1
r0p
r0d
r0a
.
Podemos assim resumir o metodo de pontos interiores primal-dual afim escala:
dado um ponto interior (x0, y0, z0) onde (x0, z0) > 0
Para k = 0, 1, · · ·
rkp = b − Axk
rkd = c − AT yk − zk
rka = −XkZke
dk = J−1(xk, yk, zk)rk
calcule o tamanho do passo primal αkp e dual αk
d
xk+1 = xk + αkpd
kx
yk+1 = yk + αkdd
ky
zk+1 = zk + αkdd
kz
6
ate convergir.
O criterio de convergencia usado neste metodo e baseado nas condicoes de otimalidade:
∣
∣
∣
xtz1+ctx+bty
∣
∣
∣ ≤ ε
‖b − Ax‖
1 + ‖b‖≤ ε,
‖c − Aty − z‖
1 + ‖c‖≤ ε,
onde ε e a tolerancia desejada e xtz uma medida da condicao de complementaridade.
O tamanho do passo αkp, αk
d e calculado da seguinte forma:
primeiramente calculamos:
ρkp = min
dxki <0
−xk
i
dxki
, (2.1)
e em seguida,
ρkd = min
dzki <0
−zk
i
dzki
. (2.2)
Os valores ρkp e ρk
d representam o tamanho de passo maximo tal que a primeira variavel
de x e z se anulam respectivamente.
Portanto, o tamanho do passo sera obtido multiplicando ρ por um parametro τ ∈
(0,1):
αkp = min(1, τρk
p) e
αkd = min(1, τρk
d)(2.3)
o que garante que nenhuma variavel de x ou z se anule. Adicionalmente, o passo e limitado
ao tamanho maximo 1, que corresponde a um passo de newton.
2.3.1 Metodo Primal-Dual Classico
O metodo afim-escala tem uma desvantagem importante: ele permite que as variaveis
(x,z) se aproximem de seus limites muito rapidamente. Consequentemente as direcoes
calculadas sao muito distorcidas e o metodo converge lentamente, pois xizi ≈ 0. Para
evitar que isto ocorra e acrescentada uma perturbacao (µ) na condicao de complemen-
taridade xizi = 0 [27]. Em seu lugar consideramos xizi = µ ou seja, o metodo primal-dual
7
resolve o seguinte sistema de equacoes nao lineares a cada iteracao:
b − Axk = 0
c − Atyk − zk = 0
µke − XkZke = 0,
onde µk e um parametro que varia a cada iteracao µk → 0 quando k → ∞. As unicas
alteracoes do metodo primal-dual classico em relacao ao metodo afim-escala sao a substi-
tuicao de rka por rk
c = µke − XkZke e o calculo da pertubacao µk = σk γk
n, onde γk = xtz,
σ ∈ (0,1), o valor γk
nseria uma “medida” da distancia media de um ponto a uma solucao
otima. E importante ressaltar que o Jacobiano permanece o mesmo o que significa que
ambos metodos necessitam resolver sistemas lineares com a mesma matriz.
2.4 Metodo Preditor-Corretor
O metodo preditor-corretor [13] e o metodo de pontos interiores mais utilizado na pratica,
pois obtem os melhores resultados computacionais e tem convergencia quadratica con-
siderando a resolucao de dois sistemas lineares em uma unica iteracao.
Este metodo calcula uma direcao com tres componentes [27]:
• direcao afim-escala,
• direcao de centragem (σ),
• direcao de correcao que compensa a aproximacao linear do metodo de Newton, que
pode ser visualizada pela relacao (x + dx)t(z + dz) = dt
xdz.
No metodo preditor corretor calcula-se primeiro a direcao afim escala(d = dx, dy, dz):
Adx = rp
Atdy + dz = rd
Zdx + Xdz = ra = −XZe
(2.4)
Usando o mesmo Jacobiano encontra-se em seguida a direcao perturbada no ponto,
(x, y, z) = (x, y, z) + (dx, dy, dz)
8
ou seja:
Adx = rp = b − Ax = b − A(x + dx) = 0
Atdy + dz = z − AT y − z = 0
Zdx + Xdz = µe − (X + Dx)(Z + Dz)e = µe − DxDze
onde, µe e a direcao de centragem e (X + Dx)(Z + Dz) representa a correcao nao linear.
A direcao utilizada sera a soma das duas direcoes:
dx = dx + dx,
dy = dy + dy,
dz = dz + dz.
Logo, (dx, dy, dz), podem ser calculados diretamente somando o dois sistemas lineares:
Adx = rp
Atdy + dz = rd
Zdx + Xdz = ra + µe − DxDze = rs.
(2.5)
ou seja, nao e necessario que d = (dx, dy, dz) seja calculada.
O calculo da perturbacao µk e funcao da direcao afim [13]. Quanto melhor a direcao
menor sera a perturbacao e vice-versa: considere,
µk = σk γk
n; γk = (xk + αk
pdkx)
t(zk + αkdd
kz)
onde, αkp e αk
d sao o tamanho do passo em relacao ao problema primal e dual, respectiva-
mente para a direcao afim, calculados da seguinte forma:
ρkp = min
dxki <0
−xi
k
dxik
e ρkd = min
dzki <0
−zi
k
dzik, sendo αk
p = min(1, τ ρkp) e αk
d = min(1, τ ρkd),
com τ ∈ (0, 1). E,
σk =
( γk
γk )3, se γk > 1
γk√
n, caso contrario.
(2.6)
A ideia desse metodo e calcular a direcao afim-escala e estudar o progresso do metodo ao
longo dessa direcao, calculando a perturbacao µk de acordo com este progresso. Uma vez
que uma segunda direcao e calculada, tambem calcula-se a correcao nao linear utilizando
o mesmo Jacobiano da direcao anterior [13]. Por outro lado, se γk < 1 existem razoes
teoricas para calcular uma perturbacao proporcional a (γk)2 [25], explicando assim a
segunda relacao no calculo de σ.
9
Podemos resumir o metodo preditor-corretor da seguinte forma:
Dados τ ∈ (0, 1), (x0, y0, z0) tal que (x0, z0) > 0
para k = 0, 1, · · ·
rkp = b − Axk
rkd = c − AT yk − zk
rka = −XkZke
Calcule dk usando (2.4)
Calcule αkp, α
kd
γk = (xk + αkpd
kx)
t(zk + αkdd
kz)
Calcule µk usando (2.6)
rks = µke + rk
a − DxkDzke
Calcule dk usando (2.5)
Calcule αkp, α
kd
xk+1 = xk + αkpd
kx
yk+1 = yk + αkdd
ky
zk+1 = zk + αkdd
kz
ate convergir.
O criterio de convergencia e o mesmo dos metodos primais-duais.
2.5 Calculo das Direcoes nos Metodos de Pontos In-
teriores
O sistema linear (2.5) tem dimensao 2n+m. Este sistema pode ser reduzido a dimensao
m atraves das seguintes relacoes:
dy = (AD−1AT )−1(rp + AD−1rd − AZ−1rs)
dx = D−1(AT dy − rd + X−1rs)
dz = X−1(rs − Zdx),
obtidas eliminando-se dz e dx sucessivamente onde D = X−1Z.
Essas simplificacoes tambem sao validas para os metodos primal-dual afim escala
e classico. A unica alteracao sera do resıduo rs que representa os valores ra e rc, respecti-
vamente. O maior esforco computacional destes metodos, portanto, esta no calculo de dy,
10
onde uma matriz simetrica definida positiva (AD−1AT ) deve ser decomposta a cada ite-
racao. Usualmente e utilizada a decomposicao de Cholesky [7] na resolucao deste sistema
linear, obtendo-se excelente desempenho computacional [1, 2, 11, 12, 15].
11
12
Capıtulo 3
Um Modelo para Formulacao do
Tratamento por Radioterapia
Neste capıtulo sao apresentados alguns modelos de programacao linear para o tratamento
por radioterapia.
O primeiro passo no desenvolvimento de um plano de tratamento e definir a lo-
calizacao do tumor e da regiao de risco (estrutura crıtica). Em seguida, a dose do raio
e calculada individualmente conforme cada paciente. Uma formulacao para programacao
linear consiste em minimizar a dose total emitida sujeito a um limite inferior da dose na
regiao do tumor. Este modelo pode ser representado pela seguinte formulacao:
min f =∑m
i=1
∑m
j=1 Dij
s. a Dij =∑n
p=1 ωpApij ∀ (i, j)
l ≤ Dij ∀ (i, j) ∈ Γ
w ≥ 0,
onde,
Dij representa a dosagem total recebida pelo paciente na posicao (i, j);
Aij representa a dosagem nominal emitida pelo raio p na posicao (i, j);
w representa a ponderacao da dosagem emitida pelos raios;
l representa a dosagem mınima a ser recebida pelo tumor;
Γ representa o subconjunto de posicoes onde se encontra o tumor;
n representa o numero de raios;
13
Figura 3.1: Uma representacao de aparelho de tratamento para radioterapia
m representa o numero total de pixels.
Definindo os seguintes subconjuntos:
• R representa o subconjunto de posicoes onde se encontra tecido crıtico (saudavel)
ou estrutura crıtica;
• N representa o subconjunto de posicoes onde se encontram os demais tecidos saudaveis;
e possıvel formular um modelo com ponderacoes em cada regiao do paciente privilegiando
ou penalizando algumas areas quanto a quantidade de dosagem a ser recebida.
Em [23] sao apresentados outros modelos de otimizacao linear que consistem em
pequenas variacoes do modelo acima, tendo em mente o aparelho mostrado na Figura 3.1.
Suponha uma matriz de pixels de imagem n × n, e que os angulos avaliados sao
θ1, θ2, θ3, · · · , θρ. Alem disso, supomos que cada angulo esta composto por η sub-raios. Os
sistemas de tratamento modernos sao capazes de realizar combinacoes complexas entre
estes sub-raios. A geometria de um modelo usando raios elementares, onde n=2, ρ = 4 e
η = 4, e indicada na Figura 3.2. Os intervalos entre os angulos e π4, e o angulo inicial e
(usualmente) zero.
Seja x(a,i) a dose ao longo do i-esimo sub-raio do a-esimo angulo (a = 1, 2, . . . , p)
e d(p,a,i) a distancia entre o sub-raio x(a,i) e a imagem do pixel p. Definimos A(p,a,i) como
o produto de e−µd(p,a,i) e a area geometrica comum a ambos os sub-raios x(a,i) e o pixel p.
Por exemplo, na Figura 3.2 o raio hachurado correspondente a x(1,2) atinge a metade do
pixel 3 e a distancia deste pixel a este raio elementar e 3√
22
(supondo que cada pixel tenha
a mesma largura igual a 1). Consequentemente, A(3,1,2) = 12e−
3√
22
µ. Os componentes da
14
1 2
3 4
x
x
xx
x
x
x
x
x
x
x
x
x
x
(2,4)
(2,3)
(2,2)
(2,1) (1,4)
(1,3)
(1,2)
x(1,1)
(4,4)
(4,3)
(4,1)x
(3,4)
(3,3)
(3,2)
(3,1)
(4,2)
Figura 3.2: Geometria de uma imagem de pixel 2 × 2 com angulosπ4, 3π
4, 5π
4e 7π
4.
matriz de propagacao dos raios, denotada por A, sao A(p,a,i), onde as linhas de A sao
indexadas por p e as colunas sao indexadas por (a,i). O fator e−µd(p,a,i) mede atenuacao
da radiacao sobre o tecido, e os valores deste coeficiente de atenuacao, (µ), dependem
da energia do raio. O coeficiente de atenuacao define a caracterıstica do tecido, ou seja,
com relacao a sua densidade. Podemos tomar como exemplo o Raio-X: quanto menor a
densidade do tecido menor o coeficiente de atenuacao, portanto queimara mais o filme e
consequentemente este ficara mais escuro; quanto maior a densidade do tecido maior o
coeficiente de atenuacao, portanto queimara menos o filme ficando mais claro.
A matriz dose de propagacao sem atenuacao (µ=0) do exemplo da Figura 3.2, e
dada por:
0 0 12
12
0 12
12
0 12
12
0 0 0 12
12
0
0 12
12
0 12
12
0 0 0 12
12
0 0 0 12
12
0 12
12
0 0 0 12
12
0 12
12
0 12
12
0 0
12
12
0 0 0 12
12
0 0 0 12
12
0 12
12
0
.
Os elementos de x podem ser ordenados da seguinte forma:
[x(1,1) x(1,2) · · · x(1,r) x(2,1) · · · x(2,r) · · · x(ρ,1) · · · x(ρ,r)], assim a dose de radiacao no
pixel p e o p-esimo componente de Ax.
Embora a matriz de propagacao da dose permita modelar facilmente os limites
superior e inferior da radiacao para algum pixel na imagem, elaborar um plano de trata-
mento nao e uma tarefa trivial. Por exemplo, um plano de tratamento pode prescrever
15
que o tumor nao receba menos do que 80 Gy, onde Gy e a dose absorvida, usualmente
medida em joules por quilograma (J/kg), tambem denominada gray (Gy): 1Gy = 1 J/kg.
Por outro lado, o plano de tratamento pode nao permitir que alguma estrutura crıtica
receba mais do que 40 Gy.
A maioria das pesquisas tem utilizado modelos de otimizacao com restricoes linea-
res, com uma das duas funcoes objetivos mais evidente, maximizacao da dose no tumor ou
minimizacao da dose na estrutura crıtica. Uma vez que esta maximizacao leva geralmente
a altas doses, outros trabalhos buscam maximizar a dose mınima do tumor ou minimizar
a dose maxima da estrutura crıtica [10].
As metas listadas abaixo indicam que este problema tem uma grande quantidade
de parametros a considerar na decisao do que seria desejavel para um plano de tratamento:
• Transmitir uma dose uniformemente letal na regiao do tumor.
• Transmitir uma radiacao tao pequena quanto possıvel na estrutura crıtica.
• Obter uma dose total tao pequena quanto possıvel.
• Reduzir a frequencia de doses altas fora da regiao do tumor.
• Controlar o numero de raios utilizados no plano de tratamento.
O objetivo mais comum adotado na pratica consiste em transmitir a maior radiacao
possıvel no tumor [23]. Ha duas razoes para evitar este objetivo. Normalmente altos
nıveis de radiacao podem conduzir uma grande soma de necrose, e o corpo humano tem
dificuldade na eliminacao de um grande volume de tecido morto. Alem disso, celulas
doentes estao distribuıdas entre tecidos saudaveis. Consequentemente, uma dose letal
uniformemente distribuıda na regiao do tumor e crucial para o sucesso do plano de trata-
mento, pois, uma dose inferior permite que a celula cancerosa sobreviva enquanto uma
dose superior pode ter efeitos altamente indesejaveis nos tecidos vizinhos.
3.1 Formulacao do Modelo
Suponha que cada orgao do corpo humano esta dividido por pixels, onde cada pixel
representara parte do tecido saudavel ou do tumor em um orgao doente. Seja, mT o
16
numero de pixels do tumor, mC o numero de pixels da estrutura crıtica e mG o numero de
pixels restantes (tecido saudavel), ou seja, m = mG +mT +mC . Finalmente n, representa
o numero de sub-raios para atingir o alvo. A prescricao e definida por quatro limites:
• ut: representa o vetor de limite superior para radiacao no tumor (ut ∈ <mT );
• lt: representa o vetor de limite inferior para radiacao no tumor (lt ∈ <mT );
• uc: representa o vetor de limite superior para radiacao na estrutura crıtica (uc ∈
<mC);
• ug: representa o vetor de limite superior para radiacao no restante de tecido saudavel
(ug ∈ <mG).
Fazemos as suposicoes obvias que 0 < lt ≤ ut, uc ≥ 0, e ug ≥ 0. Se uma dose letal
uniforme e transmitida ao tumor, o limite superior e inferior para os pixels do tumor
sao obtidos atraves das metas estabelecidas. Supondo que as metas estabelecidas para
uma celula cancerosa sejam tg, os valores de uti e lti sao geralmente (1 + ε)tg e (1 − ε)tg,
respectivamente, onde, ε e a porcentagem da variacao para a dosagem do tumor e e
denominado nıvel de uniformidade do tumor. Valores tıpicos de ε encontrados na literatura
vao de 0.02 a 0.15 [8]. O vetor ug representa a maior quantidade de radiacao permitida
para algum pixel (saudavel). Em geral tecidos saudaveis nao devem receber mais do que
10% da dose estabelecida para o tumor. Ou seja, ug = tg(1 + 0.10).
As linhas da matriz de propagacao da dose podem ser reordenadas considerando
as linhas que correspondem ao tumor, as linhas que correspondem a estrutura crıtica, e
as linhas que correspondem ao tecido saudavel. Esta reordenacao e representada pelas
sub-matrizes AT , AC e AG, como indicado abaixo:
A =
AT
AC
AG
ou seja, AT : tumor, AC : estrutura crıtica e AG : restante de tecido saudavel.
Sub-raios que nao atingem o tumor sao removidos pela eliminacao das colunas de
A que tem o vetor zero na coluna correspondente da submatriz AT . Assim, sem perda de
generalidade, consideraremos que a matriz AT nao tem colunas nulas. Portanto, temos
que A ∈ <m×n, AT ∈ <mT×n, AC ∈ <mC×n e AG ∈ <mG×n.
17
3.2 Problema de Otimizacao Linear
Um novo modelo de programacao linear usado para auxiliar no plano de radioterapia foi
introduzido em [8]. Este modelo incorpora restricoes elasticas, e quando resolvidas pelo
metodo de pontos interiores, produzem planos favoraveis. A funcao objetivo e represen-
tada pela soma ponderada de tres metas: lT t, que mede o quanto falta para que o plano
encontrado consiga aplicar a dose mınima na regiao do tumor; uTc c que mede a quantidade
de radiacao acima da prescrita recebida pela regiao crıtica; e uTg g que mede a quantidade
de radiacao acima da prescrita nos demais tecidos saudaveis. O escalar positivo w pon-
dera a importancia da formulacao de um plano que obtenha a dose mınima na regiao do
tumor, isto e, valores grandes de w forcam lT t a ser tao pequeno quanto possıvel. Seria
desejavel que existisse valor para um finito w > 0 tal que o valor otimo da componente
lT t fosse zero o que garantiria ao tumor receba o nıvel mınimo de radiacao necessario para
sua eliminacao. O modelo proposto pode ser representado pela seguinte formulacao:
min wlT t + uTc c + uT
g g
s.a lt − Lt ≤ AT x ≤ ut
ACx ≤ uc + UCc
AGx ≤ ug + UGg
0 ≤ Lt ≤ lt
−uc ≤ UCc
UGg ≥ 0
x ≥ 0.
(3.1)
onde,
w: escalar positivo,
x: dose do sub-raio que entra na imagem para alcancar o pixel p, (x ∈ <n);
t: t ∈ <mT , t ≥ 0;
c: c ∈ <mC ;
g: g ∈ <mG , g ≥ 0;
As restricoes lt−Lt ≤ AT x, ACx ≤ uc+UCc, e AGx ≤ ug+UGg, sao denominadas
elasticas, pois seus limites podem variar de acordo com os vetores t, c, e g, respectivamente.
As matrizes L, UC e UG definem como medir a elasticidade, e l, uc e ug controlam a
18
penalizacao ou recompensa com relacao a elasticidade. Valores fixos de l, uc, ug, L, UC
e UG definem um conjunto de funcoes elasticas. E estas sao incorporadas pelas seguintes
razoes: 1) a restricao elastica garante que algum conjunto de funcoes elasticas, (3.1) e
sempre estritamente factıvel; 2) a diferenca dos limites inferiores nas funcoes elasticas nos
permitem incorporar diferentes objetivos de tratamento.
Considerando que conjuntos diferentes de funcoes elasticas determinam diferentes
filosofias de tratamento a interpretacao do modelo (3.1) depende da escolha deste conjunto.
Em [8] foram propostas as seguintes escolhas, analise media e analise absoluta.
Na analise media,
l = 1mT
e, onde l ∈ <mT ;
uc = 1mC
e, onde uC ∈ <mC ;
uG = 1mG
e, onde uG ∈ <mG ; (3.2)
L = I, onde L ∈ <mT×mT ;
UC = I, onde UC ∈ <mC×mC ;
UG = I, onde UG ∈ <mG×mG .
sendo I a matriz identidade. Esta escolha tem os seguintes objetivos:
• minimizar a dosagem media recebida pelo tumor dentro do limite prescrito;
• minimizar a dosagem media da radiacao que a estrutura crıtica recebe, e
• minimizar a dosagem media que o tecido saudavel recebe.
Na analise absoluta,
l = emT, onde l ∈ <mT ;
uc = emC, onde uC ∈ <mC ;
uG = emG, onde uG ∈ <mG ; (3.3)
L = emTeT
mT, onde L ∈ <mT ;
UC = emCeT
mC, onde UC ∈ <mC ;
UG = emGeT
mG, onde UG ∈ <mG .
Esta escolha tem os seguintes objetivos:
19
• minimizar a dosagem maxima recebida pelo tumor dentro do limite prescrito;
• minimizar a dosagem maxima da radiacao que a estrutura crıtica recebe, e
• minimizar a dosagem maxima recebida pelo tecido saudavel.
3.3 Propriedades das Restricoes Elasticas
Considere a seguinte definicao:
Definicao A prescricao (ut, lt, uc, ug) admite a uniformidade do tratamento do tumor se
existe um plano, x ≥ 0, tal que lt ≤ AT x ≤ ut.
O Teorema abaixo foi demonstrado em [8].
Teorema 3.3.1 Seja (x∗(w), t∗(w), c∗(w), g∗(w)) uma solucao otima de (3.1). Para qual-
quer conjunto de funcoes elasticas temos que lT t∗(w) = O( 1w), desde que a prescricao
admita a uniformidade do tratamento do tumor.
Este Teorema mostra que se a prescricao do tratamento e viavel, o deficit de radiacao no
tumor e uniformemente limitada pelo inverso de w. Utilizando este resultado e possıvel,
resolvendo apenas um problema de programacao linear, interpretar o resultado e concluir
se existe um tratamento que atende a prescricao violando ou nao os limites ug e uc [8].
3.3.1 Interpretacao das Solucoes
Em [8] Holder mostra que quando w = κε
e o valor otimo de lT t e menor do que ε, entao
a uniformidade na regiao do tumor e garantida, onde κ e uma constante que depende
dos dados do problema. Portanto, dados w como definido acima e solucao do problema
linear correspondente, podemos interpretar a solucao para as analises media e absoluta
da seguinte forma:
Analise Media
(l = e, uc = e, ug = e)
Caso 1 : lT t∗(w) > ε, concluımos que a prescricao nao permite a uniformidade do tumor.
Caso 2 : lT t∗(w) ≤ ε, concluımos que a prescricao permite a uniformidade do tumor.
20
Esta situacao contem dois importante subcasos:
Caso 2a: uTc c∗(w) + uT
g g∗(w) > 0, a conclusao e que a uniformidade do tumor e acessıvel,
mas e cara, pois tecidos saudaveis estao recebendo mais radiacao do que desejado.
Caso 2b: uTc c∗(w)+uT
g g∗(w) ≤ 0, a conclusao e que a uniformidade do tumor e permitida,
a soma media de radiacao sobre o tecido saudavel e no mınimo tao boa quanto desejada.
Analise Absoluta
(l = uc = ug = 1)
Caso 1 : t∗(w) > ε, mesma interpretacao da analise media.
Caso 2 : t∗(w) ≤ ε, idem. Esta situacao contem dois importante subcasos:
Caso 2a: uTc c∗(w) + uT
g g∗(w) > 0, mesma interpretacao da analise media.
Caso 2b: uTc c∗(w) + uT
g g∗(w) ≤ 0, idem.
21
22
Capıtulo 4
Metodos de Pontos Interiores
Aplicados ao Modelo de Tratamento
por Radioterapia
Neste capıtulo aplicaremos os metodos de pontos interiores para o modelo (3.1) con-
siderando duas situacoes a analise media e analise absoluta.
4.1 Analise Media
Adotamos inicialmente a analise media (3.2), descrita no Capıtulo 3. Com esta escolha o
problema (3.1) se reduz a:
min w 1mT
eT t + 1mC
eT c + 1mG
eT g
s.a lt − t ≤ AT x ≤ ut
ACx ≤ uc + c
AGx ≤ ug + g
0 ≤ t ≤ lt
−uc ≤ c
(g, x) ≥ 0.
Nosso proximo objetivo consiste em obter um problema de programacao linear na forma
padrao, contendo apenas restricoes de igualdade. Seja, AT x = a, substituindo a restricao
23
canalizada por lt − t ≤ a ≤ ut e uc + c = c ⇒ c = c − uc, teremos entao:
min w 1mT
eT t + 1mC
eT (c − uc) + 1mG
eT g
s.a lt − t ≤ a ≤ ut
AT x = a
ACx ≤ c
AGx ≤ ug + g
0 ≤ t ≤ lt
(c, g, x) ≥ 0.
Introduzimos as variaveis de folga e desconsideramos a constante eT uc, obtendo o problema
primal:
min w 1mT
eT t + 1mC
eT c + 1mG
eT g
s.a a + su = ut
a + t − sl = lt
AT x − a = 0
ACx + sc − c = 0 (4.1)
AGx − g + sg = ug
t + st = lt
(t, c, g, x, su, sl, sc, sg, st) ≥ 0,
onde sem perda de generalidade substituımos c por c. Primeiramente escrevemos a matriz
de restricoes de (4.1):
M =
I 0 0 0 0 I 0 0 0 0
I 0 I 0 0 0 −I 0 0 0
−I AT 0 0 0 0 0 0 0 0
0 AC 0 −I 0 0 0 I 0 0
0 AG 0 0 −I 0 0 0 I 0
0 0 I 0 0 0 0 0 0 I
Encontraremos agora o problema dual. Multiplicando a transposta de M, por y, o vetor
de variaveis duais temos:
24
M ty =
I I −I 0 0 0
0 0 AtT At
C AtG 0
0 I 0 0 0 I
0 0 0 −I 0 0
0 0 0 0 −I 0
I 0 0 0 0 0
0 −I 0 0 0 0
0 0 0 I 0 0
0 0 0 0 I 0
0 0 0 0 0 I
yu
yl
ya
yc
yg
yt
ou seja,
yu + yl − ya = 0
AtT ya + At
Cyc + AtGyg ≤ 0
yl + yt ≤ w 1mT
e
−yc ≤1
mCe
−yg ≤ 1mG
e
yu ≤ 0
−yl ≤ 0
yc ≤ 0
yg ≤ 0
yt ≤ 0.
Definindo,
yu = −yu; yc = −yc; yg = −yg; yt = −yt e
25
substituindo no sistema anterior, obtemos:
−yu + yl − ya = 0
AtT ya − At
C yc − AtGyg ≤ 0
yl − yt ≤ w 1mT
e
yc ≤1
mCe
yg ≤ 1mG
e
(yu, yl, yc, yg, yt) ≥ 0.
Isso implica que:
−yu + yl − ya = 0
AtT ya − At
C yc − AtGyg ≤ 0
yl − yt ≤ w 1mT
e
0 ≤ yc ≤1
mCe
0 ≤ yg ≤ 1mG
e
(yu, yl, yt) ≥ 0.
Introduzindo agora as variaveis de folga z e w temos:
−yu + yl − ya = 0
AtT ya − At
C yc − AtGyg + zx = 0
yl − yt + zt = w 1mT
e
yc + wc = 1mC
e
yg + wg = 1mG
e
(yu, yl, yc, yg, yt) ≥ 0
(zx, zt) ≥ 0
(wc, wg) ≥ 0.
26
Sem perda de generalidade, tomamos yu = yu; yc = yc; yg = yg; yt = yt e, encontramos
finalmente o problema dual no formato desejado:
max −uttyu + lttyl − ut
gyg − lttyt
s.a −yu + yl − ya = 0
AtT ya − At
Cyc − AtGyg + zx = 0
yl − yt + zt = w 1mT
e (4.2)
yc + wc = 1mC
e
yg + wg = 1mG
e
(yu, yl, yc, yg, yt, zx, zt, wc, wg) ≥ 0.
As condicoes de complementaridade [27] para os problemas (4.1) e (4.2) sao dadas por:
YuSue = 0; YlSle = 0
YcSce = 0; YgSge = 0; YtSte = 0
XZxe = 0; TZte = 0;
CWce = 0; GWge = 0.
(4.3)
As condicoes de otimalidade para os problemas (4.1) e (4.2) sao dadas pela factibilidade
primal:
a + su = ut
a + t − sl = lt
AT x − a = 0
ACx + sc − c = 0
AGx − g + sg = ug
t + st = lt
(t, c, g, x, su, sl, sc, sg, st) ≥ 0,
factibilidade dual:
−yu + yl − ya = 0
AtT ya − At
Cyc − AtGyg + zx = 0
yl − yt + zt = w 1mT
e
yc + wc = 1mC
e
yg + wg = 1mG
e
(yu, yl, yc, yg, yt, zx, zt, wc, wg) ≥ 0,
e as condicoes de complementaridade (4.3).
27
4.1.1 O Metodo Primal-Dual para Analise Media
Agora, seguindo os passos descritos no Capıtulo 2, determinamos as direcoes do metodoprimal-dual, para os problemas (4.1) e (4.2) escrevendo primeiramente o Jacobiano:
J =
I 0 0 0 0 I 0 0 0 0 0 0 0 0 0 0 0 0 0 0
I 0 I 0 0 0 −I 0 0 0 0 0 0 0 0 0 0 0 0 0
−I AT 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 AC 0 −I 0 0 0 I 0 0 0 0 0 0 0 0 0 0 0 0
0 AG 0 0 −I 0 0 0 I 0 0 0 0 0 0 0 0 0 0 0
0 0 I 0 0 0 0 0 0 I 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 −I I −I 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 AtT −At
C −AtG 0 I 0 0 0
0 0 0 0 0 0 0 0 0 0 0 I 0 0 0 −I 0 I 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 I 0 0 0 0 I 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 I 0 0 0 0 I
0 0 0 0 0 Yu 0 0 0 0 Su 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 Yl 0 0 0 0 Sl 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 Yc 0 0 0 0 0 Sc 0 0 0 0 0 0
0 0 0 0 0 0 0 0 Yg 0 0 0 0 0 Sg 0 0 0 0 0
0 0 0 0 0 0 0 0 0 Yt 0 0 0 0 0 St 0 0 0 0
0 Zx 0 0 0 0 0 0 0 0 0 0 0 0 0 0 X 0 0 0
0 0 Zt 0 0 0 0 0 0 0 0 0 0 0 0 0 0 T 0 0
0 0 0 Wc 0 0 0 0 0 0 0 0 0 0 0 0 0 0 C 0
0 0 0 0 Wg 0 0 0 0 0 0 0 0 0 0 0 0 0 0 G
Multiplicando esta matriz pelas seguintes direcoes:
dt = (da, dx, dt, dc, dg, dsu, dsl
, dsc, dsg
, dst, dyu
, dyl, dya
, dyc, dyg
, dyt, dzx
, dzt, dwc
, dwg), e igua-
lando aos resıduos: r = (r1, r2, r3, · · · , r20), onde:
28
r1 = ut − a − su
r2 = lt − a − t + sl
r3 = a − AT x
r4 = −ACx − sc + c
r5 = ug − AGx + g − sg
r6 = lt − t − st
r7 = yu − yl + ya
r8 = −AtT ya + At
Cyc + AtGyg − zx
r9 = w 1mT
e − yl + yt − zt
r10 = 1mC
e − yc − wc
r11 = 1mG
e − yg − wg
r12 = µe − YuSuemT
r13 = µe − YlSlemT
r14 = µe − YcScemC
r15 = µe − YgSgemG
r16 = µe − YtStemT
r17 = µe − XZxe
r18 = µe − TZtemT
r19 = µe − CWcemC
r20 = µe − GWgemG,
(4.4)
29
obtemos assim o sistema linear, Jd=r:
da + dsu= r1 (4.5)
da + dt − dsl= r2 (4.6)
−da + AT dx = r3 (4.7)
ACdx − dc + dsc= r4 (4.8)
AGdx − dg + dsg= r5 (4.9)
dt + dst= r6 (4.10)
−dyu+ dyl
− dya= r7 (4.11)
AtT dya
− AtCdyc
− AtGdyg
+ dzx= r8 (4.12)
dyl− dyt
+ dzt= r9 (4.13)
dyc+ dwc
= r10 (4.14)
dyg+ dwg
= r11 (4.15)
Yudsu+ Sudyu
= r12 (4.16)
Yldsl+ Sldyl
= r13 (4.17)
Ycdsc+ Scdyc
= r14 (4.18)
Ygdsg+ Sgdyg
= r15 (4.19)
Ytdst+ Stdyt
= r16 (4.20)
Zxdx + Xdzx= r17 (4.21)
Ztdt + Tdzt= r18 (4.22)
Wcdc + Cdwc= r19 (4.23)
Wgdg + Gdwg= r20. (4.24)
4.1.2 Eliminacao de Variaveis
As equacoes de (4.5) a (4.24) definem o sistema linear que determina as direcoes dos
metodos de pontos interiores. Este sistema pode ser resolvido diretamente mas esta
opcao e muito cara computacionalmente devido a sua grande dimensao. No entanto, este
sistema pode ser reduzido atraves da substituicao de variaveis de forma similar a reducao
realizado para problema na forma padrao.
30
Isolando as variaveis duais de (4.16) a (4.24), obtemos:
dyu= S−1
u (r12 − Yudsu)
dyl= S−1
l (r13 − Yldsl)
dyc= S−1
c (r14 − Ycdsc)
dyg= S−1
g (r15 − Ygdsg)
dyt= S−1
t (r16 − Ytdst)
dzx= X−1(r17 − Zxdx) (4.25)
dzt= T−1(r18 − Ztdt)
dwc= C−1(r19 − Wcdc)
dwg= G−1(r20 − Wgdg)
substituindo em (4.11) a (4.15):
−S−1u (r12 − Yudsu
) + S−1l (r13 − Yldsl
) − dya= r7 ⇒
S−1u Yudsu
− S−1l Yldsl
− dya= r7 + S−1
u r12 − S−1l r13 = r21
AtT dya
− AtCS−1
c (r14 − Ycdsc) − At
GS−1g (r15 − Ygdsg
) + X−1(r17 − Zxdx) = r8 ⇒
AtT dya
+AtCS−1
c Ycdsc+At
GS−1g Ygdsg
−X−1Zxdx = r8+AtCS−1
c r14+AtGS−1
g r15−X−1r17 = r22
S−1l (r13 − Yldsl
) − S−1t (r16 − Ytdst
) + T−1(r18 − Ztdt) = r9 ⇒
− S−1l Yldsl
+ S−1t Ytdst
− T−1Ztdt = r9 − S−1l r13 + S−1
t r16 − T−1r18 = r23
S−1c (r14 − Ycdsc
) + C−1(r19 − Wcdc) = r10 ⇒
− S−1c Ycdsc
− C−1Wcdc = r10 − S−1c r14 − C−1r19 = r24
S−1g (r15 − Ygdsg
) + G−1(r20 − Wgdg) = r11 ⇒
− S−1g Ygdsg
− G−1Wgdg = r11 − S−1g r15 − G−1r20 = r25.
O novo sistema linear e dado pelas equacoes de (4.5) a (4.10) em conjunto com as equacoes
31
abaixo:
S−1u Yudsu
− S−1l Yldsl
− dya= r21 (4.26)
AtT dya
+ AtCS−1
c Ycdsc+ At
GS−1g Ygdsg
− X−1Zxdx = r22 (4.27)
−S−1l Yldsl
+ S−1t Ytdst
− T−1Ztdt = r23 (4.28)
−S−1c Ycdsc
− C−1Wcdc = r24 (4.29)
−S−1g Ygdsg
− G−1Wgdg = r25. (4.30)
Tomemos agora:
dsu= r1 − da
dsl= −r2 + da + dt
dsc= Y −1
c Sc(−C−1Wcdc − r24) (4.31)
dsg= Y −1
g Sg(−G−1Wgdg − r25)
dst= r6 − dt
e substituindo as equacoes acima em (4.8), (4.9) e (4.26) a (4.28) temos:
Acdx − dc + Y −1c Sc(−C−1Wcdc − r24) = r4 ⇒
Acdx − (I + Y −1c ScC
−1Wc)dc = r4 + Y −1c Scr24 = r26
Agdx − dg + Y −1g Sg(−G−1Wgdg − r25) = r5 ⇒
Agdx − (I + Y −1g SgG
−1Wg)dg = r5 + Y −1g Sgr25 = r27
S−1u Yu(r1 − da) − S−1
l Yl(−r2 + da + dt) − dya= r21 ⇒
− (S−1u Yu + S−1
l Yl)da − S−1l Yldt − dya
= r21 − S−1u Yur1 − S−1
l Ylr2 = r28
AtT dya
+AtCS−1
c Yc(Y−1c Sc(−C−1Wcdc−r24))+At
GS−1g Yg(Y
−1g Sg(−G−1Wgdg−r25))−X−1Zxdx =
r22 ⇒ AtT dya
− AtCC−1Wcdc − At
GG−1Wgdg − X−1Zxdx = r22 + AtCr24 + At
Gr25 = r29
− S−1l Yl(−r2 + da + dt) + S−1
t Yt(r6 − dt) − T−1Ztdt = r23 ⇒
− S−1l Ylda − (S−1
l Yl + S−1t Yt + T−1Zt)dt = r23 − S−1
l Ylr2 − S−1t Ytr6 = r30.
32
Restando as seguintes equacoes:
−da + AT dx = r3
Acdx − (I + Y −1c ScC
−1Wc)dc = r26
Agdx − (I + Y −1g SgG
−1Wg)dg = r27
−(S−1u Yu + S−1
l Yl)da − S−1l Yldt − dya
= r28
AtT dya
− AtCC−1Wcdc − At
GG−1Wgdg − X−1Zxdx = r29
−S−1l Ylda − (S−1
l Yl + S−1t Yt + T−1Zt)dt = r30.
(4.32)
Definimos agora as seguintes matrizes diagonais para simplificar a notacao:
D1 = I + Y −1c ScC
−1Wc
D2 = I + Y −1g SgG
−1Wg
D3 = S−1u Yu + S−1
l Yl
D4 = S−1l Yl (4.33)
D5 = C−1Wc
D6 = G−1Wg
D7 = X−1Zx
D8 = S−1l Yl + S−1
t Yt + T−1Zt.
Substituindo no sistema (4.32), temos:
−da + AT dx = r3 (4.34)
Acdx − D1dc = r26 (4.35)
Agdx − D2dg = r27 (4.36)
−D3da − D4dt − dya= r28 (4.37)
AtT dya
− AtCD5dc − At
GD6dg − D7dx = r29 (4.38)
−D4da − D8dt = r30. (4.39)
Isolando dyae dt das equacoes (4.37)e (4.39), teremos:
dya= −D3da − D4dt − r28 e (4.40)
dt = −D−18 (r30 + D4da). (4.41)
33
Substituımos na Equacao (4.38), obtemos:
AtT (−D3da − D4dt − r28) − At
CD5dc − AtGD6dg − D7dx = r29 ⇒
− AtT D3da − At
T D4dt − AtCD5dc − At
GD6dg − D7dx = r29 + AtT r28 ⇒
− AtT D3da − At
T D4(−D−18 (r30 + D4da)) − At
CD5dc − AtGD6dg − D7dx = r29 + At
T r28 ⇒
− AtT D3da + At
T D4D−18 D4da − At
CD5dc − AtGD6dg − D7dx = r29 + At
T r28 − AtT D4D
−18 r30
−AtT (D3 + D4D
−18 D4)da − At
CD5dc − AtGD6dg − D7dx = r31. (4.42)
Logo, temos o seguinte sistema linear:
AT dx − da = r3
ACdx − D1dc = r26
AGdx − D2dg = r27
−AtT Dada − At
CD5dc − AtGD6dg − D7dx = r31,
onde,
Da = (D3 − D4D−18 D4). (4.43)
Lembrando que:
A =
AT
AC
AG
e definindo d =
da
dc
dg
, D =
I 0 0
0 D1 0
0 0 D2
, rp =
r3
r26
r27
e
D =
Da 0 0
0 D5 0
0 0 D6
.
Podemos reescrever o sistema da seguinte forma:
Adx − Dd = rp (4.44)
−D7dx − AtDd = r31. (4.45)
Isolando d da equacao (4.44), temos:
d = D−1(Adx − rp). (4.46)
Admitindo-se que a matriz A tem mais linhas (que estao representadas pelo numero de
pixels) do que colunas (que sao representadas pelo numero de raios) e substituindo (4.46)
34
em (4.45) temos:
−D7dx − AT DD−1Adx + AT DD−1rp = r31 ⇒
(D7 + AT DD−1A)dx = r31 − AT DD−1rp = −rp.
Logo,
(AT DD−1A + D7)dx = rp. (4.47)
Este sistema e simetrico e definido positivo podendo ser resolvido pela decomposicao de
Cholesky. Sua dimensao e muito menor que o sistema original representado pelas equacoes
(4.5) a (4.24).
4.1.3 Resumo do Metodo de Pontos Interiores
O metodo de pontos interiores desenvolvido nesta secao pode ser resumido da seguinte
forma:
Seja, h= (t, c, g, x, su, sl, sc, sg, st) e j= (yu, yl, yc, yg, yt, zx, zt, wc, wg)
Dado um ponto interior (h0, j0)
Para k = 1, 2, · · ·
Calcule as matrizes diagonais D1 a D8 por (4.33) e Da por (4.43)
Calcule os resıduos:
r1 a r20 por (4.4)
r21 a r25 por (4.26) a (4.30)
r26 a r30 por (4.35) a (4.39)
r31 por (4.42)
rp por (4.44) e rp por (4.47)
Calcule d por (4.47),
Calcule as direcoes:
dyapor (4.40), dt por (4.41),
dsu, dsl
, dsc, dsg
, dstpor (4.31) e
dyu, dyc
, dyg, dyt
, dzx, dzt
, dwc, dwg
por (4.25)
Calcule ρp por (2.1) e ρd por (2.2)
Calcule αp e αd por (2.3)
ak+1 = ak + αpdka
35
yk+1a = yk
a + αddkya
hk+1p = hk + αk
pdkh
jk+1d = jk + αk
ddkj
ate convergir.
4.1.4 Metodo Preditor-Corretor para Analise Media
No metodo preditor-corretor dois sistemas lineares determinam as direcoes. Primeira-
mente e calculada a direcao afim (dh, dj, da, dya) resolvendo o sistema linear (4.5) a (4.24)
com µ = 0. Em seguida, a direcao desejada e obtida resolvendo o seguinte sistema linear
[13]:
da + dsu= r1
da + dt − dsl= r2
−da + AT dx = r3
ACdx − dc + dsc= r4
AGdx − dg + dsg= r5
dt + dst= r6
−dyu+ dyl
− dya= r7
AtT dya
− AtCdyc
− AtGdyg
+ dzx= r8
dyl− dyt
+ dzt= r9
dyc+ dwc
= r10
dyg+ dwg
= r11
Yudsu+ Sudyu
= r12
Yldsl+ Sldyl
= r13
Ycdsc+ Scdyc
= r14
Ygdsg+ Sgdyg
= r15
Ytdst+ Stdyt
= r16
Zxdx + Xdzx= r17
Ztdt + Tdzt= r18
Wcdc + Cdwc= r19
Wgdg + Gdwg= r20
36
onde os novos resıduos sao dados por:
r12 = µemT− YuSuemT
− DyuDsuemT
r13 = µemT− YlSlemT
− DylDslemT
r14 = µemC− YcScemC
− DycDscemC
r15 = µemG− YgSgemG
− DygDsgemG
r16 = µemT− YtStemT
− DytDstemT
r17 = µe − ZxXe − DzxDxe
r18 = µemT− ZtTemT
− DztDtemT
r19 = µemC− WcCemC
− DwcDcemC
r20 = µemG− WgGemG
− DwgDgemG
e o parametro µ e calculado da forma descrita no Capıtulo 2.
4.2 Analise Absoluta
Trabalhamos agora com a analise absoluta (3.3), descrita no Capıtulo 3. Com esta escolha
o problema (3.1) se reduz a:
min weTmT
t + eTmC
c + eTmG
g
s.a lt − emTeT
mTt ≤ AT x ≤ ut
ACx ≤ emC+ emC
eTmC
c
AGx ≤ emG+ emG
eTmG
g
0 ≤ emTeT
mTt ≤ lt
−emC≤ emC
eTmC
c
emGeT
mGg ≥ 0
x ≥ 0.
Nosso proximo objetivo consiste em obter um problema linear na forma padrao contendo
apenas restricoes de igualdade. Seja, AT x = a, e λ = min(lt) e como existem mT e mG
37
equacoes lineares, teremos entao:
min weTmT
t + eTmC
c + eTmG
g
s.a lt − emTeT
mTt ≤ a ≤ ut
AT x = a
ACx ≤ emC+ emC
eTmC
c
AGx ≤ emG+ emG
eTmG
g
0 ≤ eTmT
t ≤ λ
eTmC
c ≥ −1
eTmG
g ≥ 0
x ≥ 0.
Os vetores t, c e g nao tem significado fısico e neste modelo so nos interessa a soma dos
seus elementos. Assim, criamos uma variavel para caracterizar a soma dos elementos de
cada um deles: eTmT
t = τ, eTmC
c = γ, eTmG
g = β, obtendo o seguinte problema:
min wτ + γ + β
s.a lt − emTτ ≤ a ≤ ut
AT x = a
ACx ≤ emC+ emC
γ
AGx ≤ emG+ emG
β
0 ≤ τ ≤ λ
γ ≥ −1
β ≥ 0
x ≥ 0.
38
Tomando γ = γ + 1 para facilitar o desenvolvimento, e introduzindo as variaveis de folga
para obter as equacoes de igualdade e temos:
min wτ + γ − 1 + β
s.a a + su = ut
a + emTτ − sl = lt
AT x − a = 0
ACx − emCγ + sc = 0
AGx − emGβ + sg = emG
τ + σt = λ
(x, su, sl, sc, sg, σt, τ, β, γ) ≥ 0,
a : livre.
Retiramos o -1 da funcao objetivo e sem perda de generalidade substituımos γ no lugar
de γ obtemos o problema primal:
min wτ + γ + β
s.a a + su = ut
a + emTτ − sl = lt
AT x − a = 0
ACx − emCγ + sc = 0 (4.48)
AGx − emGβ + sg = emG
τ + σt = λ
(x, su, sl, sc, sg, σt, τ, β, γ) ≥ 0,
a : livre.
39
Para encontrar o problema dual, escrevemos a matriz de restricoes de (4.48):
M =
I 0 0 0 0 I 0 0 0 0
I emT0 0 0 0 −I 0 0 0
−I 0 0 0 AT 0 0 0 0 0
0 0 −emC0 AC 0 0 I 0 0
0 0 0 −emGAG 0 0 0 I 0
0 1 0 0 0 0 0 0 0 1
.
Multiplicando a transposta de M, por y, o vetor de variaveis duais temos:
M ty =
I I −I 0 0 0
0 eTmT
0 0 0 1
0 0 0 −eTmC
0 0
0 0 0 0 −eTmG
0
0 0 ATT AT
C ATG 0
I 0 0 0 0 0
0 −I 0 0 0 0
0 0 0 I 0 0
0 0 0 0 I 0
0 0 0 0 0 1
yu
yl
ya
yc
yg
πt
40
ou seja,
yu + yl − ya = 0
eTmT
yl + πt ≤ w
−eTmC
yc ≤ 1
−eTmG
yg ≤ 1
ATT ya + AT
Cyc + ATGyg ≤ 0
yu ≤ 0
−yl ≤ 0
yc ≤ 0
yg ≤ 0
πt ≤ 0.
Definindo, yu = −yu; yc = −yc; yg = −yg; πt = −πt
e substituindo no modelo anterior, obtemos:
−yu + yl − ya = 0
eTmT
yl − πt ≤ w
eTmC
yc ≤ 1
eTmG
yg ≤ 1
ATT ya − AT
C yc − ATGyg ≤ 0
(yu, yl, yc, yg, πt,≥ 0
ya : livre.
41
Introduzimos agora as variaveis de folga z, ζt, ζc e ζg obtemos as restricoes na forma de
igualdade:
−yu + yl − ya = 0
eTmT
yl − πt + ζt = w
eTmC
yc + ζc = 1
eTmG
yg + ζg = 1
ATT ya − AT
C yc − ATGyg + zx = 0
(yu, yl, yc, yg, πt, ζt, ζcζgzx) ≥ 0
ya : livre.
Sem perda de generalidade, tomamos yu = yu; yc = yc; yg = yg; πt = πt e, encontramos
finalmente o problema dual no formato desejado:
max −uTt yu + lTt yl − eT
mgyg − λπt
s.a −yu + yl − ya = 0
eTmT
yl − πt + ζt = w
eTmC
yc + ζc = 1 (4.49)
eTmG
yg + ζg = 1
ATT ya − AT
Cyc − ATGyg + zx = 0
(yu, yl, yc, yg, πt, ζt, ζc, ζg, zx) ≥ 0
ya : livre.
As condicoes de complementaridade para os problemas (4.48) e (4.49) sao dadas por:
YuSuemT= 0; YlSlemT
= 0; YcScemC= 0
YgSgemG= 0; πtσt = 0; τζt = 0
γζc = 0; βζg = 0; XZxe = 0.
(4.50)
42
As condicoes de otimalidade para os problemas (4.48) e (4.49) sao dadas pela factibilidade
primal:
a + su = ut
a + emTτ − sl = lt
AT x − a = 0
ACx − emCγ + sc = 0
AGx − emGβ + sg = emG
τ + σt = λ
(x, su, sl, sc, sg, σt, τ, β, γ) ≥ 0,
a : livre,
factibilidade dual:
−yu + yl − ya = 0
eTmT
yl − πt + ζt = w
eTmC
yc + ζc = 1
eTmG
yg + ζg = 1
ATT ya − AT
Cyc − ATGyg + zx = 0
(yu, yl, yc, yg, πt, ζt, ζc, ζg, zx) ≥ 0
ya : livre,
e as condicoes de complementaridade (4.50).
4.2.1 O Metodo Primal-Dual para Analise Absoluta
Seguindo os passos descritos no Capıtulo 2, determinamos as direcoes do metodo primal-dual, para os problemas (4.48) e (4.49) escrevendo primeiramente o Jacobiano:
J =
I 0 0 0 0 I 0 0 0 0 0 0 0 0 0 0 0 0 0 0
I 0 emT0 0 0 −I 0 0 0 0 0 0 0 0 0 0 0 0 0
−I AT 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 AC 0 −emC0 0 0 I 0 0 0 0 0 0 0 0 0 0 0 0
0 AG 0 0 −emG0 0 0 I 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 −I −I I 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 eTmt
0 0 −1 1 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 eTmc
0 0 0 1 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 eTmg
0 0 0 1 0
0 0 0 0 0 0 0 0 0 0 ATT 0 0 −AT
C −ATG 0 0 0 0 I
0 0 0 0 0 Yu 0 0 0 0 0 Su 0 0 0 0 0 0 0 0
0 0 0 0 0 0 Yl 0 0 0 0 0 Sl 0 0 0 0 0 0 0
0 0 0 0 0 0 0 Yc 0 0 0 0 0 Sc 0 0 0 0 0 0
0 0 0 0 0 0 0 0 Yg 0 0 0 0 0 Sg 0 0 0 0 0
0 0 0 0 0 0 0 0 0 πt 0 0 0 0 0 σt 0 0 0 0
0 0 ζt 0 0 0 0 0 0 0 0 0 0 0 0 0 τ 0 0 0
0 0 0 ζc 0 0 0 0 0 0 0 0 0 0 0 0 0 γ 0 0
0 0 0 0 ζg 0 0 0 0 0 0 0 0 0 0 0 0 0 β 0
0 Zx 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 X
43
Obtemos agora os resıduos resultantes da aplicacao do metodo de Newton as condicoes
de otimalidade:
r1 = ut − a − su
r2 = lt − a − emTτ + sl
r3 = −AT x + a
r4 = −ACx + emCγ − sc
r5 = emG− AGx + emG
β − sg
r6 = λ − τ − σt
r7 = yu − yl + ya
r8 = w − eTmT
yl + πt − ζt
r9 = 1 − eTmC
yc − ζc
r10 = 1 − eTmG
yg − ζg
r11 = −ATT ya + AT
Cyc + ATGyg − zx
r12 = −YuSuemT+ µemT
r13 = −YlSlemT+ µemT
r14 = −YcScemC+ µemC
r15 = −YgSgemG+ µemG
r16 = −πtσt + µ
r17 = −τζt + µ
r18 = −γζc + µ
r19 = −βζg + µ
r20 = −XZxe + µe.
(4.51)
Multiplicamos a matriz do Jacobiano pelas seguintes direcoes:
dt = (da, dx, dτ , dγ, dβ, dsu, dsl
, dsc, dsg
, dσt, dya
, dyu, dyl
, dyc, dyg
, dπt, dζt
, dζc, dζg
, dzx), e as igualamos
44
aos resıduos r = (r1, r2, r3, · · · , r20), temos, Jd=r:
da + dsu= r1 (4.52)
da + emTdτ − dsl
= r2 (4.53)
−da + AT dx = r3 (4.54)
ACdx − emCdγ + dsc
= r4 (4.55)
AGdx − emGdβ + dsg
= r5 (4.56)
dτ + dσt= r6 (4.57)
−dya− dyu
+ dyl= r7 (4.58)
eTmT
dyl− dπt
+ dζt= r8 (4.59)
eTmC
dyc+ dζc
= r9 (4.60)
eTmG
dyg+ dζg
= r10 (4.61)
ATT dya
− ATCdyc
− ATGdyg
+ dzx= r11 (4.62)
Yudsu+ Sudyu
= r12 (4.63)
Yldsl+ Sldyl
= r13 (4.64)
Ycdsc+ Scdyc
= r14 (4.65)
Ygdsg+ Sgdyg
= r15 (4.66)
πtdσt+ σtdπt
= r16 (4.67)
ζtdτ + τdζt= r17 (4.68)
ζcdγ + γdζc= r18 (4.69)
ζgdβ + βdζg= r19 (4.70)
Zxdx + Xdzx= r20. (4.71)
4.2.2 Eliminacao de Variaveis
As equacoes de (4.52) a (4.71) definem o sistema linear que determina as direcoes dos
metodos de pontos interiores para a analise media. Este sistema pode ser resolvido dire-
tamente mas esta opcao e muito cara computacionalmente devido a sua grande dimensao.
No entanto, este problema pode ser reduzido atraves da substituicao de variaveis.
45
Isolando as variaveis duais de (4.63) a (4.71), obtemos:
dyu= S−1
u (r12 − Yudsu)
dyl= S−1
l (r13 − Yldsl)
dyc= S−1
c (r14 − Ycdsc)
dyg= S−1
g (r15 − Ygdsg)
dπt= σ−1
t (r16 − πtdσt) (4.72)
dζt= τ−1(r17 − ζtdτ )
dζc= γ−1(r18 − ζcdγ)
dζg= β−1(r19 − ζgdβ)
dzx= X−1(r20 − Zxdx)
substituindo (4.72) em (4.58) a (4.62):
−dya− S−1
u (r12 − Yudsu) + S−1
l (r13 − Yldsl) = r7 ⇒
− dya+ S−1
u Yudsu− S−1
l Yldsl= r7 + S−1
u r12 − S−1l r13 = r21
eTmT
S−1l (r13 − Yldsl
) − σ−1t (r16 − πtdσt
) + τ−1(r17 − ζtdτ ) = r8 ⇒
− eTmT
S−1l Yldsl
+ σ−1t πtdσt
− τ−1ζtdτ = r8 + σ−1t r16 − τ−1r17 − eT
mTS−1
l r13 = r22
eTmC
S−1c (r14 − Ycdsc
) + γ−1(r18 − ζcdγ) = r9 ⇒
− eTmC
S−1c Ycdsc
− γ−1ζcdγ = r9 − eTmC
S−1c r14 − γ−1r18 = r23
eTmG
S−1g (r15 − Ygdsg
) + β−1(r19 − ζgdβ) = r10 ⇒
− eTmG
S−1g Ygdsg
− β−1ζgdβ = r10 − eTmG
S−1g r15 − β−1r19 = r24
AtT dya
− AtCS−1
c (r14 − Ycdsc) − At
GS−1g (r15 − Ygdsg
) + X−1(r20 − Zxdx) = r11 ⇒
AtT dya
+AtCS−1
c Ycdsc+At
GS−1g Ygdsg
−X−1Zxdx = r11+AtCS−1
c r14+AtGS−1
g r15−X−1r20 = r25.
O novo sistema linear e dado pelas equacoes de (4.52) a (4.57) em conjunto com as
46
equacoes abaixo:
−dya+ S−1
u Yudsu− S−1
l Yldsl= r21 (4.73)
−eTmT
S−1l Yldsl
+ σ−1t πtdσt
− τ−1ζtdτ = r22 (4.74)
−eTmC
S−1c Ycdsc
− γ−1ζcdγ = r23 (4.75)
−eTmG
S−1g Ygdsg
− β−1ζgdβ = r24 (4.76)
AtT dya
+ AtCS−1
c Ycdsc+ At
GS−1g Ygdsg
− X−1Zxdx = r25. (4.77)
Isolando dsuda equacao (4.52) e dσt
da equacao (4.57) obtemos:
dsu= r1 − da (4.78)
dσt= r6 − dτ (4.79)
e substituindo as equacoes acima em (4.73) e (4.74) temos:
S−1u Yu(r1 − da) − S−1
l Yldsl− dya
= r21 ⇒
− S−1u Yuda − S−1
l Yldsl− dya
= r21 − S−1u Yur1 = r26
− eTmT
S−1l Yldsl
+ σ−1t πt(r6 − dτ ) − τ−1ζtdτ = r22 ⇒
− eTmT
S−1l Yldsl
− σ−1t πtdτ − τ−1ζtdτ = r22 − σ−1
t πtr6 = r27.
Restando as seguintes equacoes:
da + emTdτ − dsl
= r2
AT dx − da = r3
ACdx − emCdγ + dsc
= r4
AGdx − emGdβ + dsg
= r5
−eTmT
S−1l Yldsl
− (σ−1t πt + τ−1ζt)dτ = r27
−eTmC
S−1c Ycdsc
− γ−1ζcdγ = r23
−eTmG
S−1g Ygdsg
− β−1ζgdβ = r24
AtT dya
+ AtCS−1
c Ycdsc+ At
GS−1g Ygdsg
− X−1Zxdx = r25
−S−1u Yuda − S−1
l Yldsl− dya
= r26.
(4.80)
47
Definindo as seguintes matrizes diagonais para simplificar a notacao:
D1 = S−1u Yu
D2 = S−1l Yl
D3 = σ−1t πt
D4 = τ−1ζt
D5 = S−1c Yc (4.81)
D6 = γ−1ζc
D7 = S−1g Yg
D8 = β−1ζg
D9 = X−1Zx
e substituindo (4.81)no modelo acima, temos:
da + emTdτ − dsl
= r2 (4.82)
AT dx − da = r3 (4.83)
ACdx − emCdγ + dsc
= r4 (4.84)
AGdx − emGdβ + dsg
= r5 (4.85)
−eTmT
D2dsl− (D3 + D4)dτ = r27 (4.86)
−eTmC
D5dsc− D6dγ = r23 (4.87)
−eTmG
D7dsg− D8dβ = r24 (4.88)
AtT dya
+ AtCD5dsc
+ AtGD7dsg
− D9dx = r25 (4.89)
−D1da − D2dsl− dya
= r26. (4.90)
Isolando dslda equacao (4.90) temos:
dsl= −D−1
2 (r26 + D1da + dya) (4.91)
Substituımos nas equacoes (4.82) e (4.86) obtendo:
da + emTdτ + D−1
2 (r26 + D1da + dya) = r2 ⇒
(I + D−12 D1)da + emT
dτ + D−12 dya
= r2 − D−12 r26 = r28
48
− eTmT
D2(−D−12 (r26 + D1da + dya
)) − (D3 + D4)dτ = r27 ⇒
eTmT
D1da + eTmT
dya− (D3 + D4)dτ = r27 − eT
mTr26 = r29.
Logo restaram as seguintes equacoes:
AT dx − da = r3 (4.92)
ACdx − emCdγ + dsc
= r4 (4.93)
AGdx − emGdβ + dsg
= r5 (4.94)
eTmT
D1da + eTmT
dya− (D3 + D4)dτ = r29 (4.95)
−eTmC
D5dsc− D6dγ = r23 (4.96)
−eTmG
D7dsg− D8dβ = r24 (4.97)
AtT dya
+ AtCD5dsc
+ AtGD7dsg
− D9dx = r25 (4.98)
(I + D−12 D1)da + emT
dτ + D−12 dya
= r28. (4.99)
Isolando da de (4.99), temos:
da = D10(r28 − emTdτ − D−1
2 dya) (4.100)
onde e D10 = (I + D−12 D1)
−1.Substituindo da em (4.92) e (4.95) teremos:
AT dx − D10(r28 − emTdτ − D−1
2 dya) = r3 ⇒
AT dx + D10emTdτ + D10D
−12 dya
= r3 + D10r28 = r30
eTmT
D1D10(r28 − emTdτ − D−1
2 dya) + eT
mTdya
− (D3 + D4)dτ = r29 ⇒
− eTmT
(D1D10D−12 − I)dya
− (eTmT
D1D10emT+ (D3 + D4))dτ = r29 − eT
mTD1D10r28 = r31.
Logo restaram as seguintes equacoes:
AT dx + D10emTdτ + D10D
−12 dya
= r30 (4.101)
ACdx − emCdγ + dsc
= r4 (4.102)
AGdx − emGdβ + dsg
= r5 (4.103)
−eTmT
(D1D10D−12 − I)dya
− (eTmT
D1D10emT+ (D3 + D4))dτ = r31 (4.104)
−eTmC
D5dsc− D6dγ = r23 (4.105)
−eTmG
D7dsg− D8dβ = r24 (4.106)
AtT dya
+ AtCD5dsc
+ AtGD7dsg
− D9dx = r25 (4.107)
49
Lembrando que,
A =
AT
AC
AG
e definindo
d1 =
dτ
dγ
dβ
, d2 =
dya
dsc
dsg
, E =
D10emT0 0
0 −emC0
0 0 −emG
,
D =
D10D−12 0 0
0 I 0
0 0 I
,
D =
−(eTmT
D1D10emT+ (D3 + D4)) 0 0
0 −D6 0
0 0 −D8
,
E =
−eTmT
(D1D10D−12 − I) 0 0
0 −eTmC
D5 0
0 0 −eTmG
D7
,
rp =
r30
r4
r5
, rp =
r31
r23
r24
, D =
I 0 0
0 D5 0
0 0 D7
.
Assim substituımos (4.101), (4.102) e (4.103) por:
Adx + Ed1 + Dd2 = rp.
Agora restaram as seguintes equacoes:
Adx + Ed1 + Dd2 = rp (4.108)
Ed2 + Dd1 = rp (4.109)
AtDd2 − D9dx = r25. (4.110)
50
Isolando d1 de (4.109) teremos:
d1 = D−1(rp − Ed2). (4.111)
Substituindo em (4.108) temos:
Adx + ED−1(rp − Ed2) + Dd2 = rp ⇒
Adx − ED−1Ed2 + Dd2 = rp − ED−1rp = rs.
Assim obtemos:
Adx − (ED−1E − D)d2 = rs (4.112)
AtDd2 − D9dx = r25. (4.113)
Isolando agora d2 de (4.112) temos:
d2 = (ED−1E − D)−1(Adx − rs) (4.114)
Substituindo em (4.113) teremos:
AtD(ED−1E − D)−1(Adx − rs) − D9dx = r25 ⇒
AtD(ED−1E − D)−1Adx − D9dx = r25 + AtD(ED−1E − D)−1rs = r.
Logo restou a seguinte equacao:
(AtD(ED−1E − D)−1A + D9)dx = r. (4.115)
Este sistema e simetrico e definido positivo podendo ser resolvido pela decomposicao de
Cholesky. Sua dimensao e muito menor que o sistema original representado pelas equacoes
(4.52) a (4.71).
A primeira vista, a inversao da matriz (D−ED−1E)−1 parece ser computacional-
mente cara. Mas se utilizarmos a formula de Sherman-Morrison-Woodbury [6] obtemos a
seguinte relacao:
(D − ED−1E)−1 = D−1 + D−1ED−1(I3 − ED−1ED−1)−1ED−1 (4.116)
que envolve a inversao de uma matriz diagonal I3 − ED−1ED−1 de dimensao 3. Portanto
esta inversa pode ser facilmente calculada nao representando grande esforco computa-
cional.
51
4.2.3 Resumo do Metodo de Pontos Interiores
O metodo de pontos interiores desenvolvido nesta secao pode ser resumido da seguinte
forma:
Seja, f= (x, su, sl, sc, sg, σt, τ, β, γ) e g = (yu, yl, yc, yg, πt, ζt, ζc, ζg, zx)
Dado um ponto interior (f 0, g0)
Para k = 1, 2, · · ·
Calcule as matrizes diagonais D1 a D9 por (4.81)
Calcule os resıduos:
r1 a r20 por (4.51)
r21 a r25 por (4.73) a (4.77)
r26 e r27 por (4.90)e(4.86)
r28 e r29 por (4.99) e (4.95)
r30 e r31 por (4.101) e (4.104)
rp e rp por (4.108) e (4.109)
rs e r por (4.112) e (4.115)
Calcule d2 por (4.114),
Calcule d1 por (4.111),
Calcule as direcoes:
da por (4.100), dslpor (4.91),
dsupor (4.78), dσt
por (4.79)
dyu, dyl
, dyc, dyg
, dπt, dζt
, dζc, dζg
, dzxpor (4.72)
Calcule αp e αd de forma equivalente a descrita no Capıtulo 2
ak+1 = ak + αpdka
yk+1a = yk
a + αpdkya
fk+1 = fk + αkpd
kf
gk+1 = gk + αkdd
kg
ate convergir.
52
4.2.4 Metodo Preditor-Corretor para Analise Absoluta
No metodo preditor-corretor dois sistemas lineares determinam as direcoes. Primeira-
mente e calculada a direcao afim (df , dg, da, dya) resolvendo o sistema linear (4.52) a
(4.71) com µ = 0. Em seguida, a direcao desejada e obtida resolvendo o seguinte sistema
linear [13]:
da + dsu= r1
da + emTdτ − dsl
= r2
−da + AT dx = r3
ACdx − emCdγ + dsc
= r4
AGdx − emGdβ + dsg
= r5
dτ + dσt= r6
−dya− dyu
+ dyl= r7
eTmT
dyl− dπt
+ dζt= r8
eTmC
dyc+ dζc
= r9
eTmG
dyg+ dζg
= r10
ATT dya
− ATCdyc
− ATGdyg
+ dzx= r11
Yudsu+ Sudyu
= r12
Yldsl+ Sldyl
= r13
Ycdsc+ Scdyc
= r14
Ygdsg+ Sgdyg
= r15
πtdσt+ σtdπt
= r16
ζtdτ + τdζt= r17
ζcdγ + γdζc= r18
ζgdβ + βdζg= r19
Zxdx + Xdzx= r20
53
onde os novos resıduos sao dados por
r12 = µemT− YuSuemT
− DyuDsuemT
r13 = µemT− YlSlemT
− DylDslemT
r14 = µemC− YcScemC
− DycDscemC
r15 = µemG− YgSgemG
− DygDsgemG
r16 = µ − πtσt − dπtdσt
r17 = µ − ζtτ − dζtdτ
r18 = µ − ζcγ − dζcdγ
r19 = µ − ζgβ − dζgdβ
r20 = µe − ZxXe − DzxDxe
e o parametro µ e calculado da forma descrita no Capıtulo 2.
54
Capıtulo 5
Resultados Computacionais
Neste capıtulo as versoes de metodos de pontos interiores primal-dual e preditor-corretor
para as analises media e absoluta sao comparadas. Os metodos foram implementados
em Matlab 5.3 e os testes realizados em um microcomputador PC compatıvel 900MHz
usando o sistema operacional Linux.
Nos experimentos foram utilizados dois problemas. O primeiro e um exemplo teste
pequeno de apenas 16 pixels e o segundo com 4096 pixels baseados em um problema real
obtido na pagina www.trinity.edu/aholder/research/oncology/.
As Tabelas 5.1 a 5.4 resumem os resultados obtidos quanto ao numero de iteracoes,
contagem de operacoes de ponto flutuante (flops ou megaflops) e tempo de execucao em
segundos respectivamente, para os metodos primal-dual e preditor-corretor utilizando
tanto o comando interno do Matlab quanto a decomposicao de Cholesky. A tolerancia
utilizada para convergencia e a raiz quadrada do ε da maquina [5] e o parametro τ tem o
valor fixo 0,99995.
Na analise absoluta foram implementadas duas versoes. A primeira delas utiliza
a matriz do lado esquerdo da equacao (4.116), a segunda versao (rapida) utiliza o lado
direito da mesma equacao. Somente a versao rapida foi utilizada nos experimentos para
o problema de 4096 pixels.
Na Tabela 5.1 podemos ver que ambos os metodos se comportam bem para um
problema pequeno e como seria de se esperar o metodo primal dual converge mais rapido
pois o preditor corretor foi desenvolvido para problemas de grande porte. Da mesma
forma, a utilizacao da decomposicao de Cholesky no matlab para este problema nao
55
proporciona ganhos significativos.
Metodo Iteracoes Flops Tempo
Primal-Dual com Cholesky 12 81099 0,05
Primal-Dual sem Cholesky 12 84075 0,05
Preditor-Corretor com Cholesky 10 91568 0,07
Preditor-Corretor sem Cholesky 10 98548 0,06
Tabela 5.1: Exemplo com 16 pixels - Analise Media
Na Tabela 5.2 os resultados para a analise absoluta sao similares ao da analise
media, sendo que a utilizacao da decomposicao de Cholesky obteve convergencia ligeira-
mente mais rapida. A versao rapida que inverte a matriz de forma eficiente reduziu
quase a metade o numero de flops necessario a convergencia e tambem reduz o numero
de iteracoes.
Metodo Iteracoes Flops Tempo
Primal-Dual com Cholesky 9 123034 0,06
Primal-Dual sem Cholesky 9 125266 0,06
Preditor-Corretor com Cholesky 6 109103 0,05
Preditor-Corretor sem Cholesky 7 131980 0,06
Primal-Dual rapido com Cholesky 8 65713 0,05
Primal-Dual rapido sem Cholesky 8 67697 0,05
Preditor-Corretor rapido com Cholesky 6 61445 0,05
Preditor-Corretor rapido sem Cholesky 6 65645 0,05
Tabela 5.2: Exemplo com 16 pixels - Analise Absoluta
Para o problema de grande porte (Tabela 5.3) na analise media, podemos ver que o
metodo preditor-corretor e o uso da decomposicao de Cholesky ja apresentam resultados
significativamente melhores. O matlab mascara um pouco este resultado no tempo com-
putacional mas no numero de flops podemos observar melhor este desempenho superior.
O metodo primal dual nao converge para este problema na analise absoluta (Tabela
5.4). O preditor corretor no entanto converge em um numero muito pequeno de iteracoes,
embora com um tempo computacional elevado. Estes experimentos foram realizados em
56
Metodo Iteracoes MFlops Tempo
Primal-Dual com Cholesky 14 69,0 75,3
Primal-Dual sem Cholesky 14 93,3 76,8
Preditor-Corretor com Cholesky 13 73,1 71,7
Preditor-Corretor sem Cholesky 13 141,1 75,6
Tabela 5.3: Exemplo com 4096 pixels - Analise Media
uma estacao Sun Blade 100 com o Matlab 6.0, assim, nao foi possıvel estimar o numero
de flops.
Metodo Iteracoes Tempo
Primal-Dual rapido com Cholesky Max.
Primal-Dual rapido sem Cholesky Max.
Preditor-Corretor rapido com Cholesky 8 974,3
Preditor-Corretor rapido sem Cholesky 8 975,2
Tabela 5.4: Exemplo com 4096 pixels - Analise Absoluta
57
58
Capıtulo 6
Conclusoes
Os resultados de nossa implementacao em Matlab indicam que os metodos de pontos
interiores sao promissores para esta classe de problemas sendo o metodo preditor-corretor
mais eficiente, principalmente no aspecto robustez. Pode-se observar que as iteracoes
do metodo sao rapidas devido a reducao do sistema linear via eliminacao de variaveis,
resultando em um sistema da ordem do numero de linhas, que determina o numero de
pixels, da matriz de propagacao da radiacao, obtendo em ambos os casos um sistema de
muito menor dimensao, simetrico e definido positivo que pode ser decomposto utilizando
Cholesky acelerando a convergencia do metodo.
Na analise absoluta houve uma simplificacao do modelo para facilitar o desen-
volvimento do metodo e da resolucao dos sistemas lineares. Foram criadas variaveis que
representam a soma dos elementos dos vetores cujos valores individuais nao sao impor-
tantes.
Alem disso, e possıvel reduzir o esforco computacional, utilizando a formula de
Sherman-Morrison-Woodbury para inverter uma matriz. Devido a utilizacao desta formula,
somente matrizes diagonais sao invertidas.
6.1 Perspectivas Futuras
Pretende-se comparar a implementacao desenvolvida com a implementacao descrita em
[8] para o metodo de pontos interiores, que nao realiza a reducao do sistema linear. Nesta
comparacao serao utilizados problemas reais e espera-se uma significativa reducao no
59
tempo computacional da versao especializada dos metodos.
Para obter melhor desempenho computacional explorando totalmente as estruturas
matriciais envolvidas, deve-se realizar uma implementacao utilizando linguagens como C
ou Fortran.
60
Bibliografia
[1] Ilan Adler, Narendra Karmarkar, Mauricio G. C. Resende, and Geraldo Veiga. Data
structures and programming techniques for the implementation of Karmarkar’s algo-
rithm. ORSA Journal on Computing, 1:84–106, 1989.
[2] Ilan Adler, Mauricio G. C. Resende, Geraldo Veiga, and Narendra Karmarkar. An
implementation of Karmarkar’s algorithm for linear programming. Mathematical
Programming, 44:297–335, 1989.
[3] G. K. Bahr, J. G. Kereiakes, H. Horwitz, R. Finney, J. Galvin, and K. Good. The
method of linear programming applied to radiation treatment planning. Radiology,
91:686–693, 1968.
[4] Jordi Castro. A specialized interior-point algorithm for multcommodity network
flows. SIAM J. Optimization, 10(3):852–877, 2000.
[5] John E. Dennis and Robert B. Schnabel. Numerical Methods for Unconstrained
Optimization and Nonlinear Equations. SIAM, Philadelphia, PA, 1996.
[6] I. S. Duff, A. M. Erisman, and J. K. Reid. Direct Methods for Sparse Matrices.
Clarendon Press, Oxford, 1986.
[7] Gene H. Golub and Charles F. Van Loan. Matrix Computations. Johns Hopkins - 3a
edicao, 1996.
[8] Allen Holder. Designing radiotherapy plans with elastic constraints and interior point
methods. A ser publicado em Health Care and Management Science, 6(1), 2003.
61
[9] M. Langer, R. Brown, M. Urie, J. Leong, M. Stracher, and J. Shapiro. Large scale
optimization of beam weights under dose-volume restrictions. International J. of
Radiation Oncology, Biology, Physics, 18:887–893, 1990.
[10] W. Lodwick, S. McCout, F. Newman, and S. Humphries. Series in Applied Mathe-
matics - Computational, Radiology and imaging: Therapy and Diagnosis. Optimiza-
tion methods for radiation therapy plans. In C. Borges and F. Natterer, editors, IMA.
Springer-Verlag., 1998.
[11] I. J. Lustig, R. E. Marsten, and D. F. Shanno. On implementing Mehrotra’s predictor-
corrector interior point method for linear programming. SIAM J. Optimization,
2:435–449, 1992.
[12] K. A. McShane, C. L. Monma, and D. F. Shanno. An implementation of a primal-
dual interior-point method for linear programming. ORSA Journal on Computing,
1:70–83, 1989.
[13] S. Mehrotra. On the implementation of a primal-dual interior point method. SIAM
Journal on Optimization, 2(4):575–601, 1992.
[14] Renato D. C. Monteiro, Ilan Adler, and Mauricio G. C. Resende. A polynomial-time
primal-dual affine scaling algorithm for linear and convex quadratic programming
and its power series extension. Mathematics of Operations Research, 15:191–214,
1990.
[15] Aurelio Ribeiro Leite Oliveira and Christiano Lyra. Interior point methods for the
polynomial L∞ fitting problem. Trabalho submetido ao Internacional Transactions
in Operational Research, 2002.
[16] Aurelio Ribeiro Leite Oliveira, Mario A. Nascimento, and Christiano Lyra. Efficient
implementation and benchmark of interior point methods for the polynomial L1
fitting problem. Statistics & Data Analysis, 35(2):119–135, 2000.
[17] Aurelio Ribeiro Leite Oliveira, Leonardo Nepomuceno, and Secundino Soares. Short
term hydroelectric scheduling combining network flow and interior point approaches.
Trabalho submetido a Electrical Power & Energy Systems, 2001.
62
[18] Aurelio Ribeiro Leite Oliveira, Leonardo Nepomuceno, and Secundino Soares. Op-
timal active power dispatch combining network flow and interior point approaches.
Aceito para publicacao, IEEE Transactions on Power Systems, 2003.
[19] Aurelio Ribeiro Leite Oliveira and Secundino Soares. Metodos de pontos interiores
para problema de fluxo de potencia otimo. Anais do XIII Congresso Brasileiro de
Automatica, em CD-ROM, Florianopolis, SC, pages 790–795, 2000.
[20] Aurelio Ribeiro Leite Oliveira and Secundino Soares. Metodos de pontos interiores
para problema de fluxo de potencia otimo. Trabalho aceito para publicacao na Revista
SBA: Controle & Automacao, 2003.
[21] Mauricio G. C. Resende and Geraldo Veiga. An efficient implementation of a network
interior point method. DIMACS Series in Discrete Mathematics and Theoretical
Computer Scienc e, 12:299–348, 1993.
[22] I. Rosen, R. Lane, S. Morril, and J. Belli. Treatment plan optimization using linear
programming. Med. Phys., 18:141–152, 1990.
[23] D. Shepard, M. Ferris, G. Olivera, and T. Mackie. Optimizing the delivery of radia-
tion therapy to cancer patients. SIAM Review, 41(4):721–744, 1999.
[24] D. Sonderman and P. Abrahamson. Radiotherapy treatment design using mathe-
matical programming models. Operations Research, 33(4):705–725, 1985.
[25] R. A. Tapia and Yin Zhang. Superlinear and quadratic convergence of primal-dual
interior point methods for linear programming revisited. Journal of Optimization
Theory and Applications, 73:229–242, 1992.
[26] R.J. Vanderbei. An interior-point code for quadratic programming. Technical Report
SOR-94-15, Department of civil Engineering and Operations Research, Princeton
University, Princeton, N.J., 1994.
[27] S. J. Wright. Primal–Dual Interior–Point Methods. SIAM Publications, SIAM,
Philadelphia, PA, USA, 1996.
63