Planejamento do tratamento por radioterapia através de métodos de pontos interiores

Planejamento do tratamento por radioterapia através de métodos de pontos interiores

Cecília Bollini Barboza^{*, I}; Aurelio Ribeiro Leite de Oliveira^II

^IInstituto de Ciências Matemáticas e de Computação Universidade de São Paulo (USP) São Carlos – SP [e-mail]

^IIInst. de Matemática, Estatística e Computação Científica Universidade Estadual de Campinas (UNICAMP) Campinas – SP [e-mail]

RESUMO

Neste trabalho são desenvolvidos o estudo e implementação de métodos de pontos interiores específicos para o problema de planejamento do tratamento de câncer por radioterapia. Este é um problema de grande porte que contém uma estrutura matricial particular. Esta estrutura é explorada de forma eficiente obtendo um sistema linear de dimensão muito menor. Além disso, o sistema obtido é esparso, simétrico e definido positivo. Resultados numéricos em Matlab mostram que problemas de grande porte podem ser resolvidos em poucas iterações e baixo tempo computacional por esta abordagem.

Palavras-chave: métodos de pontos interiores; programação linear; radioterapia.

ABSTRACT

In this work, a specialized interior point method is developed for planning cancer treatment 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. Moreover, the system obtained is sparse, symmetric and positive definite. Matlab numerical results show that large-scale problems can be solved in few iterations and short computational time by this approach.

Keywords: interior point methods; linear programming; radiotherapy.

1. Introdução

O tratamento por radioterapia vem sendo utilizado desde há muito tempo em pacientes que sofrem diversas formas de câncer. O objetivo é eliminar as células do tumor através de radiação ao mesmo tempo em que procura evitar a destruição de células vizinhas saudáveis. Mais da metade de pacientes com câncer são submetidos à radiação em algum ponto durante o curso de tratamento da enfermidade. Existem atualmente máquinas sofisticadas inspiradas na tomografia computadorizada que emitem radiação ao longo de todo o corpo do paciente [15], ampliando as possibilidades de minimização dos efeitos colaterais do tratamento uma vez que as diferentes alternativas de atingir o tumor podem evitar a exposição de regiões críticas à radiação e distribuir doses menores ao redor do corpo focadas na região do tumor. Do ponto de vista matemático, o desafio consiste em emitir uma alta dosagem de radiação no tumor, suficiente para sua eliminação, e simultaneamente, minimizar a radiação nas regiões vizinhas, compostas de tecido saudável, reduzindo ao máximo as complicações nestas regiões que são muitas vezes críticas.

Uma vantagem deste problema, do ponto de vista computacional, é que ele pode ser modelado como um problema de otimização linear [1,4,5,13,14,15]. Nas soluções básicas (vértices), em que as variáveis não-básicas estão em seus limites, parte do tumor recebe a dosagem mínima estabelecida para sua eliminação enquanto parte do tecido saudável pode receber a dosagem máxima permitida o que não corresponde ao tratamento mais adequado do ponto de vista médico. Na presença de múltiplas soluções os métodos de pontos interiores convergem para aquelas em certa medida distanciadas de vértices. Estas seriam as soluções mais adequadas na formulação de um tratamento, pois menos variáveis tendem a atingir seus limites na presença de múltiplas soluções. Uma motivação para a abordagem por métodos de pontos interiores está na estrutura matricial bem definida deste problema. Experiências anteriores com o desenvolvimento de métodos de pontos interiores específicos para uma classe de problemas e a exploração da estrutura matricial dos sistemas lineares resultantes mostram que esta abordagem é muito bem sucedida em termos da obtenção de melhor desempenho computacional [2,8,9,10,11,12].

Este artigo está organizado da seguinte forma. A Seção 2 descreve uma formulação do tratamento por radioterapia. Na Seção 3 discutimos o problema de otimização linear para auxiliar no tratamento por radioterapia. Na Seção 4 apresentamos os métodos de pontos interiores aplicados ao modelo de tratamento por radioterapia, considerando a análise média e a análise absoluta. Os resultados computacionais são apresentados na Seção 5 e as conclusões obtidas estão resumidas na Seção 6.

2. Formulação Matemática

O primeiro passo no desenvolvimento de um plano de tratamento é definir a localização das células cancerosas (tumor) e da região de risco (estrutura crítica). Em seguida, a intensidade (dose) da radiação é calculada individualmente conforme cada paciente. Uma abordagem por programação linear consiste em minimizar a dose total emitida sujeito a um limite inferior da dose na região do tumor [14]. Ainda em [14], são apresentadas pequenas variações do modelo de programação linear tendo em mente o aparelho mostrado na Figura 1.

É possível formular um modelo com ponderações em cada região do paciente privilegiando ou penalizando algumas áreas quanto à quantidade de dosagem a ser recebida.

Suponha para exemplificação uma matriz de imagem bidimensional m x n onde cada elemento da matriz denominado pixel representa uma (pequena) parte do corpo humano, e que os ângulos avaliados são q₁, q₂,..., q_t. Além disso, considere que cada ângulo está composto por h sub-raios. Os sistemas de tratamento modernos são capazes de realizar combinações complexas entre estes sub-raios. A geometria de um modelo usando raios elementares, com m=n=2 (4 pixels), 4 ângulos t=4 e cada ângulo é composto por 4 sub-raios, é indicada na Figura 2. O intervalo entre os ângulos é , e o ângulo inicial é .

Seja x_j a dose ao longo do i-ésimo sub-raio de ângulo a com j = (a-1)h + i " i = 1,...,h; a = 1,...,t e considere d_(p,j) a distância entre o sub-raio x_j e pixel de imagem p. Definimos A_(p,,j) como o produto de e a área geométrica comum a ambos o sub-raio x_j e o pixel p. O fator mede atenuação da radiação sobre o tecido, e os valores deste coeficiente de atenuação (µ) dependem da energia do raio. O coeficiente de atenuação define a característica do tecido com relação a sua densidade. Podemos tomar como exemplo de atenuação uma chapa Raio-X de um tecido humano. Quanto menor a densidade do tecido menor o coeficiente de atenuação, portanto queimará mais o filme e conseqüentemente este ficará mais escuro; quanto maior a densidade do tecido maior o coeficiente de atenuação, portanto queimará menos o filme ficando mais claro, em geral representando tecido ósseo.

Os elementos A_(p,j) formam a matriz de propagação, denotada por A, cujas linhas são indexadas por p e colunas são indexadas por j (a,i). Assim, a dose de radiação no pixel p é dada pelo p-ésimo componente de Ax.

Por exemplo, na Figura 2 o raio pontilhado correspondente a x₂ atinge a metade do pixel 3 e a distância deste pixel a este sub-raio é (supondo que cada pixel tenha a mesma largura). Conseqüentemente, A_{(3,2) =} . Considerando o coeficiente de atenuação nulo (µ=0), quando não existe perda de energia na propagação do raio, a radiação total sobre o pixel 3 da Figura 2 é dada pela relação:

Logo, a matriz de propagação sem atenuação do exemplo é dada por:

Embora a matriz de propagação permita modelar facilmente os limites superiores e inferiores da radiação para algum pixel de imagem, elaborar um plano de tratamento não é uma tarefa trivial. Por exemplo, um plano de tratamento pode prescrever que o tumor não receba menos do que 80 Gy, onde Gy é a dose absorvida, usualmente medida em Joules por quilograma (J/kg), também denominada Gray (Gy): 1Gy = 1 J/kg. Por outro lado, o plano de tratamento pode não permitir que alguma estrutura crítica receba mais do que 40 Gy.

A maioria das pesquisas tem utilizado modelos de otimização com restrições lineares, com uma das duas funções objetivos mais evidente, maximização da dose no tumor ou minimização da dose na estrutura crítica. Uma vez que esta maximização leva geralmente a altas doses, outros trabalhos buscam maximizar a dose mínima do tumor ou minimizar a dose máxima da estrutura crítica [6].

As metas listadas abaixo indicam que este problema tem uma grande quantidade de parâmetros a considerar na decisão do que seria desejável para um plano de tratamento:

• Transmitir uma dose uniformemente letal na região do tumor.

• Transmitir uma radiação tão pequena quanto possível na estrutura crítica.

• Obter uma dose total tão pequena quanto possível.

• Reduzir a freqüência de doses altas fora da região do tumor.

• Controlar o número de raios utilizados no plano de tratamento.

O objetivo mais comum adotado na prática consiste em transmitir a maior radiação possível no tumor [14]. Existem duas razões para evitar este objetivo: células doentes estão geralmente distribuídas entre tecidos saudáveis e o corpo humano tem dificuldade na eliminação de um volume grande de tecido morto como o produzido por altas doses de radiação. Portanto, uma dose letal uniformemente distribuída na região do tumor é crucial para o sucesso do plano de tratamento, pois, uma dose inferior permite que a célula cancerosa sobreviva enquanto uma dose superior pode ter efeitos altamente indesejáveis nos tecidos vizinhos.

Vamos considerar, que uma região do corpo humano seja representado por uma rede de pixels, onde cada pixel representará parte do tecido saudável ou do tumor em um órgão doente. Seja, m_t o número de pixels do tumor, m_c o número de pixels da estrutura crítica e m_g o número de pixels restantes (tecido saudável), m = m_g + m_t + m_c. Finalmente n, representa o número de sub-raios que atinge o volume do alvo. A prescrição é definida por quatro limites:

• u_t: representa o vetor de limite superior para radiação no tumor

  ;

• l_t: representa o vetor de limite inferior para radiação no tumor

  ;

• u_c: representa o vetor de limite superior para radiação na estrutura crítica

  ;

• u_g: representa o vetor de limite superior para radiação no restante de tecido saudável

  .

Fazemos as suposições óbvias que e Se uma dose letal uniforme é transmitida ao tumor, o limite superior e inferior para os pixels do tumor são obtidos através das metas estabelecidas. Supondo que as metas estabelecidas para uma célula cancerosa sejam t_g. Os valores de são geralmente e respectivamente, onde, e é a porcentagem da variação para a dosagem do tumor e é denominado nível de uniformidade do tumor. Valores típicos de e encontrados na literatura vão de 0.02 à 0.15 [4]. O vetor u_g representa a maior quantidade de radiação permitida para algum pixel (saudável). Em geral tecidos saudáveis não devem receber mais do que 10% da dose estabelecida para o tumor. Ou seja,

As linhas da matriz de propagação podem ser reordenadas considerando as linhas que correspondem ao tumor, as linhas que correspondem à estrutura crítica, e as linhas que correspondem ao tecido saudável. Esta reordenação é representada pelas sub-matrizes A_t, A_c e A_g, como indicado abaixo:

ou seja, A_t : tumor, A_c : estrutura crítica e A_g : restante de tecido saudável.

Sub-raios que não atingem o tumor são removidos pela eliminação das colunas de A que tem o vetor zero na coluna A_t correspondente. Assim, sem perda de generalidade, consideraremos que a matriz A_t não tem colunas nulas. Portanto, temos que e

3. Problema de Otimização Linear

Nesta seção apresentamos o modelo de programação linear usado para auxiliar no plano de radioterapia introduzido em [4]. Este modelo incorpora as denominadas restrições elásticas que satisfazem todas as restrições impostas para o tratamento quando uma solução existe ou apresentam a melhor solução fora das especificações possível de acordo com uma ponderação da função objetivo. O modelo proposto pode ser representado pela seguinte formulação [4]:

onde

x: representa o vetor de dose dos sub-raios ;

t:

c:

g:

A função objetivo é representada pela soma ponderada de três metas: l^tt, que mede o quanto falta para que o plano encontrado consiga aplicar a dose mínima na região do tumor; que mede a quantidade de radiação acima da prescrita recebida pela região crítica; e que mede a quantidade de radiação acima da prescrita nos demais tecidos saudáveis.

As restrições são denominadas elásticas, pois seus limites podem variar de acordo com os vetores t, c, e g, respectivamente. As matrizes L, U_c e U_g definem como medir a elasticidade, e l, u_c e u_g controlam a penalização ou recompensa com relação à elasticidade. Valores fixos de l, u_c, u_g, L, U_c e U_g definem um conjunto de funções elásticas. As restrições elásticas garantem que para algum conjunto de funções elásticas (1) é sempre estritamente factível. Desta forma as variáveis t, c e g se assemelham ao clássico conceito de variáveis artificiais. Adicionalmente, a diferença dos limites inferiores nas funções elásticas permite incorporar diferentes objetivos de tratamento.

O escalar positivo w pondera a importância da formulação de um plano que obtenha a dose mínima na região do tumor, isto é, valores grandes de w forçam l^t t a ser tão pequeno quanto possível priorizando a aplicação da dose mínima no tumor em relação à violação da dosagem máxima permitida nos tecidos saudáveis quando não é possível satisfazer todas as especificações. Seria desejável que existisse um valor finito para w>0 tal que o valor ótimo da componente l^t t fosse zero o que garantiria ao tumor receber o nível mínimo de radiação necessário para sua eliminação.

Considerando que conjuntos diferentes de funções elásticas determinam diferentes filosofias de tratamento a interpretação do modelo (1) depende da escolha deste conjunto. Em [4] foram propostas as seguintes escolhas, análise média e análise absoluta.

Na análise média,

sendo I a matriz identidade. Esta escolha tem os seguintes objetivos:

• minimizar a dosagem média recebida pelo tumor dentro do limite prescrito;

• minimizar a dosagem média da radiação que a estrutura crítica recebe;

• minimizar a dosagem média que o tecido saudável recebe.

Na análise absoluta,

Esta escolha tem os seguintes objetivos:

• minimizar a dosagem máxima recebida pelo tumor dentro do limite prescrito;

• minimizar a dosagem máxima da radiação que a estrutura crítica recebe;

• minimizar a dosagem máxima recebida pelo tecido saudável.

3.1 Propriedades das Restrições Elásticas

Considere a seguinte definição:

DefiniçãoA prescrição (u_t, l_t, u_c, u_g) admite a uniformidade do tratamento do tumor se existe um plano, x> 0, tal que l_t<A_t x <u_t.

O Teorema abaixo foi demonstrado em [4].

Teorema 3.1 Seja (x^*(w), t^*(w), c^*(w), g^*(w)) uma solução ótima de (1). Para qualquer conjunto de funções elásticas temos que l^tt^*(w) = O, desde que a prescrição admita a uniformidade do tratamento do tumor.

Este Teorema mostra que se a prescrição do tratamento é viável, o déficit de radiação no tumor é uniformemente limitada pela inversa de w. Utilizando este resultado é possível, resolvendo apenas um problema de programação linear, interpretar o resultado e concluir se existe um tratamento que atende à prescrição violando ou não os limites u_g e u_c [4].

3.1.1 Interpretação das Soluções

Em [4] Holder mostra que quando w = e o valor ótimo de l^t t é menor do que e, então a uniformidade na região do tumor é garantida, onde k é uma constante que depende dos dados do problema. Portanto, dados w como definido acima e solução do problema linear correspondente, podemos interpretar a solução para as análises média e absoluta da seguinte forma:

Análise Média (l = u_c = u_g = e) e Análise Absoluta (l = u_c = u_g = 1)

Caso 1:l^t t^*(w) > e, concluímos que a prescrição não permite a uniformidade do tumor.

Caso 2:l^t t^*(w) < e, concluímos que a prescrição permite a uniformidade do tumor. Esta situação contém dois importantes subcasos:

Caso 2a: a conclusão é que a uniformidade do tumor é acessível, mas é cara, pois tecidos saudáveis estão recebendo mais radiação do que desejado.

Caso 2b: a conclusão é que a uniformidade do tumor é permitida, a soma de radiação sobre o tecido saudável é no mínimo tão boa quanto desejada.

4. Métodos de Pontos Interiores aplicados ao Modelo de Tratamento por Radioterapia

4.1 Análise Média

Adotamos inicialmente a análise média, descrita na Seção 3. Com esta escolha o problema (1) se reduz a:

Introduzimos as variáveis de folga e desconsideramos a constante e^t u_c, obtendo o problema primal na forma padrão:

Encontramos agora o problema dual no formato desejado:

As condições de complementaridade [16] para os problemas (4) e (5) são dadas por:

As condições de otimalidade para os problemas (4) e (5) são dadas pela factibilidade primal, factibilidade dual e as condições de complementaridade (6).

4.2 O Método Primal-Dual para Análise Média

Aplicando o método de Newton às condições de otimalidade, e modificando o tamanho do passo para manter os pontos interiores, determinamos as direções do método primal-dual, para os problemas (4) e (5). Escrevendo primeiro o Jacobiano:

Obtemos agora os resíduos resultantes da aplicação do método de Newton às condições de otimalidade:

Multiplicando a matriz do Jacobiano pelas seguintes direções:

e igualando aos resíduos r = (r₁, r₂, r₃, ..., r₂₀), temos Jd = r onde:

4.3 Eliminação de Variáveis

As equações (7) definem o sistema linear que determina as direções dos métodos de pontos interiores. Este sistema pode ser resolvido diretamente, mas esta opção é muito cara computacionalmente devido à sua grande dimensão. No entanto, este sistema pode ser reduzido através da substituição de variáveis de forma similar à redução realizado para problema na forma padrão. Admitindo-se que a matriz A tem mais linhas (que estão representadas pelo número de pixels) do que colunas (que são representadas pelo número de raios). Logo,

Este sistema é simétrico e definido positivo podendo ser resolvido pela decomposição de Cholesky. Sua dimensão é muito menor que o sistema original representado pelas equações (7). As direções são:

onde

com,

O método de pontos interiores desenvolvido nesta seção pode ser resumido da seguinte forma:

Para calcular o tamanho do passo , obtemos = - e = - .

Os valores representam o tamanho de passo máximo tal que a primeira componente de h e j se anulam respectivamente. Portanto, o tamanho do passo será obtido multiplicando por um parâmetro assim o que garante que nenhuma componente de h ou j se anule. Adicionalmente, o passo é limitado ao tamanho máximo 1, que corresponde a um passo de Newton.

4.4 Método Preditor-Corretor para Análise Média

No método preditor-corretor dois sistemas lineares determinam as direções. Primeiramente é calculada a direção afim resolvendo o sistema linear (7) com µ = 0. Em seguida, os novos resíduos são calculados:

e a direção desejada é obtida resolvendo o seguinte sistema linear [7]:

O cálculo da perturbação é função da direção afim [7]. Quanto melhor a direção menor será a perturbação e vice-versa: e onde, e são o tamanho do passo em relação à direção afim. A idéia desse método é calcular a direção afim-escala e estudar o progresso do método ao longo dessa direção, calculando a perturbação de acordo com este progresso. Uma vez que uma segunda direção é calculada, também calcula-se a correção não linear utilizando o mesmo Jacobiano da direção anterior [7].

4.5 Análise Absoluta

Trabalhamos agora com a análise absoluta (3). Com esta escolha o problema (1) se reduz a:

Os vetores t, c e g não tem significado físico e neste modelo só nos interessa a soma dos seus elementos. Assim, criamos uma variável para caracterizar a soma dos elementos de cada um deles: , obtendo o seguinte problema:

Nosso próximo objetivo consiste em obter um problema linear na forma padrão contendo apenas restrições de igualdade.

Encontramos o problema dual no formato desejado:

As condições de complementaridade para os problemas (8) e (9) são dadas por:

As condições de otimalidade para os problemas (8) e (9) são dadas pela factibilidade primal, factibilidade dual e as condições de complementaridade (10).

4.6 O Método Primal-Dual para Análise Absoluta

Aplicando o método de Newton às condições de otimalidade, e modificando o tamanho do passo para manter os pontos interiores, determinamos as direções do método primal-dual, para os problemas (8) e (9). Escrevendo primeiro o Jacobiano:

Obtemos agora os resíduos resultantes da aplicação do método de Newton às condições de otimalidade:

Multiplicamos a matriz do Jacobiano pelas seguintes direções:

e as igualamos aos resíduos temos, Jd=r:

4.7 Eliminação de Variáveis

As equações de (11) definem o sistema linear que determina as direções dos métodos de pontos interiores Para a análise média este sistema pode ser resolvido diretamente, mas esta opção é muito cara computacionalmente devido à sua grande dimensão. No entanto, este problema pode ser reduzido através da substituição de variáveis. Logo restou a seguinte equação:

Este sistema é simétrico e definido positivo podendo ser resolvido pela decomposição de Cholesky. Sua dimensão é muito menor que o sistema original representado pelas equações (11).

A primeira vista, a inversão da matriz parece ser computacionalmente cara. Mas se utilizarmos a fórmula de Sherman-Morrison-Woodbury [3] obtemos a seguinte relação:

que envolve a inversão de uma matriz diagonal de dimensão 3. Portanto esta inversa pode ser facilmente calculada não representando grande esforço computacional.

As direções são:

onde,

Também temos que,

com,

O método de pontos interiores desenvolvido nesta seção pode ser resumido da seguinte forma:

4.8 Método Preditor-Corretor para Análise Absoluta

No método preditor-corretor dois sistemas lineares determinam as direções. Primeiramente é calculada a direção afim resolvendo o sistema linear (11) com µ = 0. Em seguida, a direção desejada é obtida resolvendo o seguinte sistema linear [7]:

onde os novos resíduos são dados por

5. Resultados Computacionais

Neste capítulo as versões de métodos de pontos interiores prima-dual e preditor-corretor para as análises média e absoluta são comparadas. Os métodos 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 é um exemplo teste pequeno de apenas 16 pixels e o segundo com 4096 pixels baseados em um problema real obtido na página www.trinity.edu/aholder/research/oncology/.

As Tabelas 5.1 a 5.4 resumem os resultados obtidos quanto ao número de iterações, contagem de operações de ponto flutuante (flops ou megaflops) e tempo de execução em segundos respectivamente, para os métodos primal-dual e preditor-corretor utilizando tanto o comando interno do Matlab quanto a decomposição de Cholesky. A tolerância utilizada para convergência é a raiz quadrada do Î da máquina e o parâmetro t tem o valor fixo 0,99995.

Na análise absoluta foram implementadas duas versões. A primeira delas utiliza a matriz do lado esquerdo da equação (4.7), a segunda versão (rápida) utiliza o lado direito da mesma equação. Somente a versão rápida foi utilizada nos experimentos para o problema de 4096 pixels. Obtivemos também os seguintes resultados para as variáveis w, l^tt, u^tc, u^tg, g, c e g para a análise média e w, t, g e b para análise absoluta.

Para Análise Média:

- Primal-dual com Cholesky:

- Primal-dual sem Cholesky:

- Preditor-corretor com Cholesky:

- Preditor-corretor sem Cholesky:

Para Análise Absoluta:

w = 0.1;

t = 0.8501;

- Primal-dual com Cholesky:

g = 6.205×10^-8;

b = 5×10^-11.

- Primal-dual sem Cholesky:

g = 4.528×10^-9;

b =2.352×10^-12.

- Preditor-corretor com Cholesky:

g = 1.125×10^-7;

b = 3.81×10^-9.

- Preditor-corretor sem Cholesky:

g = 4.737×10^-10;

b = 2.2167×10^-11.

Na Tabela 5.1 podemos ver que ambos os métodos se comportam bem para um problema pequeno e como seria de se esperar o método primal dual converge mais rápido, pois o preditor corretor foi desenvolvido para problemas de grande porte. Da mesma forma, a utilização da decomposição de Cholesky no Matlab para este problema não proporciona ganhos significativos.

Na Tabela 5.2 os resultados para a análise absoluta são similares ao da análise média, sendo que a utilização da decomposição de Cholesky obteve convergência ligeiramente mais rápida. A versão rápida que inverte a matriz de forma eficiente reduziu quase à metade o número de flops necessário à convergência e também reduz o número de iterações.

Para o problema de grande porte (Tabela 5.3) na análise média, podemos ver que o método preditor-corretor e o uso da decomposição de Cholesky já apresentam resultados significativamente melhores. O Matlab mascara um pouco este resultado no tempo computacional mas no número de flops podemos observar melhor este desempenho superior.

O método primal dual não converge para este problema na análise absoluta (Tabela 5.4). O preditor corretor, no entanto converge em um número muito pequeno de iterações, embora com um tempo computacional elevado. Estes experimentos foram realizados em uma estação Sun Blade 100 com o Matlab 6.0, assim, não foi possível estimar o número de flops.

6. Conclusões

Os resultados de nossa implementação em Matlab indicam que os métodos de pontos interiores são promissores para esta classe de problemas sendo o método preditor-corretor mais eficiente, principalmente no aspecto robustez. Pode-se observar que as iterações do método são rápidas, devido à redução do sistema linear via eliminação de variáveis, resultando em um sistema da ordem do número de linhas, que determina o número de pixels, da matriz de propagação, obtendo em ambos os casos um sistema de muito menor dimensão, simétrico e definido positivo que pode ser decomposto utilizando Cholesky acelerando a convergência do método.

Na análise absoluta houve uma simplificação do modelo para facilitar o desenvolvimento do método e da resolução dos sistemas lineares. Foram criadas variáveis que representam a soma dos elementos dos vetores cujos valores individuais não são importantes.

Além disso, é possível reduzir o esforço computacional, utilizando a fórmula de Sherman- Morrison-Woodbury para inverter uma matriz. Devido à utilização desta fórmula, somente matrizes diagonais são invertidas.

AGRADECIMENTOS

Este trabalho foi parcialmente financiado pela FAPESP – Fundação de Amparo à Pesquisa do Estado de São Paulo e pelo CNPq – Conselho Nacional de Desenvolvimento Científico e Tecnológico.

Recebido em 07/2003; aceito em 10/2005 após 2 revisões

Received July 2003; accepted October 2005 after 2 revisions

• [1] Bahr, G.K.; Kereiakes, J.G.; Horwitz, H.; Finney, R.; Galvin, J. & Good, K. (1968). The Method of Linear Programming Applied to Radiation Treatment Planning. Radiology, 91, 686-693.
• [2] Castro, J. (2000). A Specialized Interior-Point Algorithm for Multcommodity Network Flows. SIAM J. Optimization, 10(3), 852-877.
• [3] Duff, I.S.; Erisman, A.M. & Reid, J.K. (1986). Direct Methods for Sparse Matrices. Clarendon Press, Oxford.
• [4] Holder, A. (2003). Designing Radiotherapy Plans with Elastic Constraints and Interior Point Methods. Health Care and Management Science, 6(1), 5-16.
• [5] Langer, M.; Brown, R.; Urie, M.; Leong, J.; Stracher, M. & Shapiro, J. (1990). Large Scale Optimization of Beam Weights under Dose-Volume Restrictions. International J. of Radiation Oncology, Biology, Physics, 18, 887-893.
• [6] Lodwick, W.; McCout, S.; Newman, F. & Humphries, S. (1998). Series in Applied Mathematics Computational, Radiology and Imaging: Therapy and Diagnosis. In: Optimization methods for radiation therapy plans [edited by C. Borges and F. Natterer], IMA, Springer-Verlag.
• [7] Mehrotra, S. (1992). On the Implementation of a Primal-Dual Interior Point Method. SIAM Journal on Optimization, 2(4), 575-601.
• [8] Oliveira, A.R.L.; Nascimento, M.A. & Lyra, C. (2000). Efficient Implementation and Benchmark of Interior Point Methods for the Polynomial L1 Fitting Problem. Computational Statistics & Data Analysis, 35, 119-135.
• [9] Oliveira, A.R.L.; Nepomuceno, L. & Soares, S. (2005). Short Term Hydroelectric Scheduling Combining Network Flow and Interior Point Approaches. Electrical Power & Energy Systems, 27, 91-99.
• [10] Oliveira, A.R.L.; Nepomuceno, L. & Soares, S. (2003). Optimal Active Power Dispatch Combining Network Flow and Interior Point Approaches. IEEE Transactions on Power Systems, 18(4), 1235-1240.
• [11] Oliveira, A.R.L. & Soares, S. (2003). Métodos de Pontos Interiores para Problema de Fluxo de Potência Ótimo DC. SBA: Controle & Automação, 14(3), 278-285.
• [12] Resende, M.G.C. & Veiga, G. (1993). An Efficient Implementation of a Network Interior Point Method. DIMACS Series in Discrete Mathematics and Theoretical Computer Science, 12, 299-348.
• [13] Rosen, I.; Lane, R.; Morril, S. & Belli, J. (1990). Treatment Plan Optimization Using Linear Programming. Med. Phys., 18, 141-152.
• [14] Shepard, D.; Ferris, M.; Olivera, G. & Mackie, T. (1999). Optimizing the Delivery of Radiation Therapy to Cancer Patients. SIAM Review, 41(4), 721-744.
• [15] Sonderman, D. & Abrahamson, P. (1985). Radiotherapy Treatment Design Using Mathematical Programming Models. Operations Research, 33(4), 705-725.
• [16] Wright, S.J. (1996). PrimalDual InteriorPoint Methods. SIAM Publications, Philadelphia PA.

*

Corresponding author / autor para quem as correspondências devem ser encaminhadas

Datas de Publicação

• Publicação nesta coleção
  01 Jun 2006
• Data do Fascículo
  Abr 2006

Histórico

• Recebido
  Jul 2003
• Aceito
  Out 2005