Desenvolvimento de um Algoritmo de Otimização Auto-Adaptativo para a Determinação de um Protocolo Otimizado para a Administração de Drogas no Tratamento de Tumores

RESUMO

Tradicionalmente, os parâmetros utilizados nos algoritmos de otimização heurísticos são considerados fixos durante o processo evolutivo. Apesar desta característica simplificar os códigos computacionais e dos bons resultados apresentados na literatura, o uso de parâmetros fixos não previne a ocorrência de convergência prematura, além de problemas relacionados à sensibilidade destes parâmetros. Neste sentido, este trabalho tem por objetivo propor um algoritmo heurístico auto-adaptativo baseado no conceito de taxa de convergência e de diversidade da população aplicados ao algoritmo de Evolução Diferencial. A metodologia proposta é aplicada na resolução de funções matemáticas com diferentes graus de complexidade e na determinação de um protocolo otimizado para a administração de drogas em pacientes com câncer através da proposição e resolução de um problema de controle ótimo multi-objetivo. Assim, deseja-se minimizar o número de células tumorais e a concentração máxima das drogas administradas ao paciente. A partir da obtenção da Curva de Pareto, um protocolo otimizado para a administração de drogas pode ser testado na prática.

Palavras-chave:
algoritmo auto-adaptativos; evolução diferencial; tratamento de tumores; problema de controle ótimo multi-objetivo

ABSTRACT

Traditionally, the parameters used in heuristic optimization algorithms areconsidered constant during the evolutionary process. Although this characteristic simplifies the computational codes and despite the good quality of results presented in the literature, the use of constant parameters does not avoid the occurrence of premature convergence and other difficulties related to parameters sensitivity. In this context, this study aims at developing a self-adaptive heuristic algorithm based on rate of convergence and population diversity concepts, which are used by the Differential Evolution algorithm. The methodology proposed is applied to the minimization of mathematical functions and to the determination of a protocol for drug administration in patients with cancer, through the formulation and solution of a multi-objective optimal control problem. In the present study the minimization of the number of cancerous cells and the minimization of the concentration of drugs that are administered to the patient represent the considered objective functions. The Pareto's Curve provides a set of optimized protocols, among which an efficient solution for drug administration can be chosen through a given criterion, aiming at practical applications.

Keywords:
self-adaptive algorithm; differential evolution algorithm; tumor treatment; multi-objective optimal control problem

1 INTRODUÇÃO

Em quaisquer métodos heurísticos de otimização é necessária a definição de parâmetros associados às diferentes técnicas para que estes possam ser aplicados. Em comparação com estratégias que fazem uso de informações sobre o gradiente da função objetivo e de suas restrições, o número de parâmetros a serem definidos pelo usuário é bem superior, o que caracteriza uma aparente desvantagem. Neste caso, existe a necessidade de se definir valores ou faixas para tais parâmetros de modo que estes possam ser utilizados para uma gama representativa de estudos de caso. Além disso, o uso de um conjunto de parâmetros default sugerido pela literatura especializada não implica na obtenção da solução global ou em um número otimizado de avaliações da função objetivo.

Tradicionalmente, parâmetros como, por exemplo, o tamanho da população, são considerados fixos durante todo o processo de otimização. Do ponto de vista de implementação, tal característica simplifica os algoritmos, mas por outro lado deixa de considerar um aspecto muito interessante da evolução biológica, a saber, a variação do número de indivíduos da população, que é inerente aos sistemas biológicos reais. Em termos práticos, se o número de gerações em um algoritmo evolutivo tender ao infinito, implica que a população caminha para se tornar homogênea, isto é, todos os candidatos da população terão o mesmo valor de função objetivo. Neste caso, se essa condição é satisfeita, não existe a necessidade de avaliar todos os candidatos, o que implica na redução do número de avaliações da função objetivo.

Neste cenário, para reduzir a pressão no processo de seleção que pode ser causada pelo uso de parâmetros constantes, nos últimos anos várias metodologias têm sido propostas para a atualização desses parâmetros. Zaharie ¹ propôs uma regra baseada no conceito de diversidade da população para a atualização dinâmica da taxa de perturbação e da probabilidade de cruzamento no algoritmo de Evolução Diferencial (ED) ^2),(3. Neste contexto, uma linha de pesquisa quetem sido bastante explorada consiste do uso de modelos caóticos de busca para a atualização destes parâmetros, devido à diversidade que esses modelos introduzem no processo evolutivo ^4),(5),(6. Em se tratanto do tamanho da população, não são muitos os trabalhos que abordam este tema. Dentre estes, Sun e colaboradores ⁷ propuseram duas estratégias para a atualização dinâmica do tamanho da população no algoritmo de Enxame de Partículas. A principal desvantagem desta estratégia é o fato destas não fazerem uso de nenhuma informação sobre o andamento do processo evolutivo, ou seja, dependem exclusivamente de fatores externos ao processo (parâmetros adotados para as constantes das funções matemáticas adotadas para a redução do número de indivíduos da população).

Em se tratando de aplicações considerando métodos de otimização heurísticos, a literatura especializada apresenta uma infinidade de estudos de caso. Dentre estes, sua utilização para o desenvolvimento de estratégias empregadas para o tratamento de tumores representa um tópico de pesquisa de grande relevância nos dias atuais. Apesar dos inúmeros esforços dedicados ao desenvolvimento de novas tecnologias para o diagnóstico e o tratamento de tumores, a incidência de câncer tem aumentado significativamente no mundo. Mesmo com as várias modalidades de tratamento, não existe um procedimento considerado ideal, isto é, aquele que é capaz de destruir as células tumorais sem prejudicar as células normais, além de inibir a manifestação da agressividade da doença ⁸. Nos dias atuais, esta busca pode ser facilitada pelo desenvolvimento de modelos matemáticos que são capazes de representar adequadamente o fenômeno de crescimento de tumores. Assim, de posse de um modelo representativo, mesmo que qualitativamente, pode-se formular um problema de otimização cujos objetivos primordiais sejam o de determinar a melhor forma de administração de drogas ao paciente, bem como simular mudanças em diferentes aspectos como a dose por fração, a taxa de dose e o tempo total do tratamento, além de identificar protocolos para a administração de drogas, na redução dos efeitos colaterais à radiação, dentre outros aspectos ^9),(10.

Matematicamente, a determinação do perfil para a administração de drogas pode ser obtida através da formulação e resolução de um Problema de Controle Ótimo (PCO). Cabe ressaltar que, tradicionalmente, o objetivo considerado neste problema é a minimização do número de células tumorais ao final do tratamento. Neste caso, não é considerado um aspecto muito relevante nestes sistemas, a saber, que altas concentrações de drogas afetam significativamente a saúde dos pacientes (efeitos colaterais evidenciais, por exemplo, a partir da redução do número de células imunes decorrente da aplicação das grandes quantidades de drogas).

Neste contexto, o objetivo deste trabalho é propor um algoritmo auto-adaptativo para a atualização dos parâmetros do algoritmo de ED através da definição do conceito de taxa de convergência e do uso do conceito de diversidade da população proposto por Zaharie ¹. A metodologia proposta será utilizada para a resolução de uma série de funções matemáticas, bem como na resolução de um PCO multi-objetivo que consiste na determinação de um conjunto de protocolos para a administração de drogas em pacientes com câncer, dentre os quais pode-se, a partir da definição de um determinado critério, implementá-lo na prática. Os objetivos considerados neste problema multi-objetivo são: a minimização do número de células tumorais e a minimização da concentração de drogas que devem ser administradas ao paciente, de modo a minimizar os efeitos colaterais que hoje não são considerados na formulação deste tipo de problema.

Este trabalho está estruturado conforme segue: na Seção 2 é apresentada a fundamentação teórica do algoritmo de ED. Já nas Seções 3 e 4 são apresentadas as motivações para o desenvolvimento deste trabalho e as estratégias para a atualização dos parâmetros considerados no algoritmo de ED. Na Seção 5 são apresentados a modelagem matemática relacionada ao crescimento de tumores, bem como a formulação matemática do PCO. Os resultados e a discussão são apresentados na Seção 6. Finalmente, as conclusões e sugestões para trabalhos futuros são apresentadas na última seção.

2 ALGORITMO DE EVOLUÇÃO DIFERENCIAL

O algoritmo de ED, proposto por Storn & Price ², é uma estratégia heurística desenvolvida para a resolução de problemas de otimização que se fundamenta na realização de operações vetoriais para a geração de candidatos em potencial para a solução do problema. O seu procedimento geral é dado pelas seguintes etapas ^2),(3:

1. Define-se o problema (número de variáveis de projeto, bem como o seu domínio, além da função objetivo e das restrições) e os parâmetros que serão utilizados no algoritmo de ED (tamanho da população - NP, taxa de perturbação - F, probabilidade de cruzamento - CR, estratégia utilizada para gerar uma nova população (ver a Tabela 1), número de gerações e critério de parada);

2. Gera-se uma população inicial com NP soluções factíveis (as variáveis de projeto estão dentro das fronteiras delimitadas pelo domínio especificado pelo usuário) para o problema em questão;

3. Seleciona-se três indivíduos (x^j _α,, x^j _β e x^j _γ , onde os subscritos α, β e γ representam diferentes posições dos indivíduos na população e j é um contador pertencente ao intervalo [1 NP]), de forma aleatória, na população (conforme a estratégia 7, por exemplo). Um deles é escolhido para ser substituído (x^j _α - candidato a ser perturbado) e os outros dois (x^j _β e x^j _γ ) para realizarem a operação de substração vetorial. Cada indivíduo é um vetor que contém todas as variáveis de projeto;

4. Gera-se uma nova população de candidatos através da obtenção de novos indivíduos (x^{j +1}, j = 1,...,NP) via realização de uma operação vetorial simples, isto é, x^{j +1} = x^j _α + F(x^j _β - x^j _γ ) (a forma com a qual esse novo candidato será gerado é função daestratégia escolhida, conforme apresentado na Tabela 1), ponderada por uma taxa de perturbação F (número real definido entre [0 2] ²). O faixa considerada para este parâmetro foi definida a partir da resolução de diferentes estudos de caso encontrados na literatura especializada ^2),(11;

5. Avalia-se essa nova população com relação à função objetivo e as restrições;

6. Compara-se a posição relativa de cada indivíduo das duas populações (original e a gerada através do mecanismo apresentado) em relação ao valor da função objetivo e se for gerado um número aleatório (pertencente ao intervalo [0 1]) e este for menor que a probabilidade de cruzamento CR definida, o novo candidato é aceito, caso contrário, esse candidato é rejeitado. Ao final deste procedimento, uma nova população (formada com indivíduos da população corrente e da nova população) será formada. Cabe ressaltar que a faixa para o parâmetro CR, isto é, [0 1] é definida por se tratar de uma probabilidade. Na prática, quanto maior esse valor, isto é, mais próximo da unidade, maior é a chance de um candidato que foi gerado substituir o seu antecessor na população;

7. Avalia-se o critério de parada definido. Se este for satisfeito o procedimento evolutivo é finalizado. Caso contrário, voltar ao terceiro passo.

Tabela 1:
Representação das estratégias do algoritmo de ED ^2),(3.

Estratégia	Representação	Mecanismo de Mutação
1	xj+1 = xj best + F(xj β - xj γ)	DE/best/1/exp
2	xj+1 = xj α + F(xj β - xj γ)	DE/rand/1/exp
3	xj+1 = xj old + F(xj best - xj old) + F(xj γ - xi δ)	DE/ran-to-best/2/exp
4	xj+1 = xj best + F(xj α - xj β) + F(xj γ - xi δ)	DE/best/2/exp
5	xj+1 = xj α + F(xj ρ - xj β) + F(xj γ - xi δ)	DE/rand/2/exp
6	xj+1 = xj best + F(xj β - xj γ)	DE/best/1/bin
7	xj+1 = xj α + F(xj β - xj γ)	DE/rand/1/bin
8	xj+1 = xj old + F(xj best - xj old) + F(xj γ - xi δ)	DE/rand-to-best/2/bin
9	xj+1 = xj best + F(xj α - xj β) + F(xj γ - xi δ)	DE/best/2/bin
10	xj+1 = xj α + F(xj ρ - xj β) + F(xj γ - xi δ)	DE/rand/2/bin

A convenção utilizada por Storn & Price ² é DE/x/y/z, onde x representa o vetor que seráperturbado. Por exemplo, optando-se por rand, a posição que contém o vetor que será perturbado é escolhido aleatoriamente na população. y é o número de pares de vetores que são considerados durante a perturbação e z é o tipo de cruzamento usado para a geração do candidato, bin no caso binomial e exp no caso exponencial. Os subscritos α, β, γ, ρ e δ são índices escolhidos aleatoriamente na população. x_best é o melhor indivíduo da população na geração anterior e x_old é um indivíduo, escolhido aleatoriamente da população, na geração anterior.

Para a inicialização do algoritmo de ED, Storn & Price ² aconselham NP em torno de 5 a 10 vezes o número de variáveis de projeto do problema, F como sendo um valor entre 0,4 a 1,2 e CR como sendo igual a 0,8. Cabe ressaltar que esses valores ou faixas podem ser diferentes a depender da aplicação considerada.

3 INFLUÊNCIA DOS PARÂMETROS DO ALGORITMO DE EVOLUÇÃO DIFERENCIAL

Conforme comentado anteriormente, a escolha dos parâmetros empregados no algoritmo de ED podem influenciar o processo de busca do ótimo global, bem como o número total de avaliações da função objetivo requeridas para a sua obtenção. Para avaliar o efeito da escolha destes parâmetros, bem como da homogeneidade da população durante o processo evolutivo, será considerada a função matemática definida pela equação (3.1) ¹².

(3.1)

em que x_i , i=1,2 é o vetor de variáveis de projeto definido no intervalo 0 ≤ x ₁, x ₂ ≤ 10 .Este problema apresenta, dentro do domínio especificado, inúmeros ótimos locais e um ótimo global cujo valor reportado na literatura é f=18,55 ¹².

Para a resolução deste problema são consideradas as seguintes configurações no algoritmo de ED: ED₁ (CR= 0,1 e F= 0,3), ED₂ (CR= 0,5 e F= 0,5) e ED₃ (CR= 0,8 e F= 1,2), onde CR é a probabilidade de cruzamento e F é a taxa de perturbação. Em cada uma destas estratégias foram considerados 50 indivíduos na população, estratégia 7 (ED/rand/1/bin - convenção utilizada em ² em que ED representa o algoritmo base, rand significa que a posição que contém o vetor a ser perturbado será escolhida aleatoriamente na população, 1 representa o número de pares que são considerados durante a perturbação e bin representa o tipo de cruzamento empregado para a geração do novo candidato). O critério de parada é definido como sendo a homogeneidade da população, isto é, o módulo da diferença entre o valor médio e o pior valor em termos da função objetivo (da geração corrente) deve ser menor que 1E-10 para que o algoritmo seja finalizado. Além disso, cada uma destas estratégias é executada dez vezes considerando um conjunto diferente de sementes ([0 1 2 ... 9]) para a inicialização do gerador de números aleatórios.

A Figura 1 apresenta a evolução do valor da função objetivo (f) dada pela equação (3.1) (CR= 0,8, F= 1,2 e semente igual a 9). Já a Figura 2 apresenta a influência da escolha do conjunto de parâmetros no número de gerações requeridas para a obtenção do ótimo global, em função da semente considerada para a inicialização do gerador de números aleatórios.

Figura 1:
Evolução da população ao longo do processo evolutivo.

Figura 2:
Influência dos parâmetros no número de gerações requeridas.

Com os resultados apresentados na Figura 1 observa-se que a partir da 50^α geração cada um dos indivíduos da população apresentam o mesmo valor para a função objetivo, isto é, a homogeneidade da população foi alcançada. Neste caso, como os candidatos da população têm o mesmo valor para a função objetivo, não existe a necessidade de se avaliar todos estes candidatos, e nem mesmo de se trabalhar com esse número de indivíduos. Intuitivamente, pode ser benéfico para a população se expandir nas primeiras gerações, isto é, assumir o seu valor máximo (50 indivíduos) enquanto a diversidade da população é alta. Este aspecto oferece aos indivíduos da população a oportunidade de exploração do espaço de projeto, promovendo assim uma melhor chance de localização do ótimo global. Por outro lado, do ponto de vista da otimização em si, ao final do processo evolutivo, como a tendência natural da população é se tornar homogênea, resultam avaliações desnecessárias da função objetivo e, por consequência, no aumento do custo computacional. Neste último caso, o tamanho da população poderá assumir um valor mínimo, o que ainda pode permitir a exploração de novas regiões, bem como o refinamento da solução corrente.

Na Figura 2 observa-se que a escolha do conjunto de parâmetros CR e F influencia o número de gerações requeridas para a obtenção do ótimo global. Neste caso, um conjunto de valores default, que é interessante para um estudo de caso, pode não ser o melhor para outras aplicações. De forma prática, a escolha destes parâmetros representa uma pressão sobre o usuário da técnica, sendo necessária a realização de inúmeros experimentos numéricos para a definição do melhor conjunto de parâmetros que conduzam à solução global, bem como leve a um número menor de avaliações da função objetivo. Além disso, conforme discutido por Zaharie ¹, apesar do desempenho e das várias aplicações considerando parâmetros fixos no algoritmo de ED, não existe nenhuma garantia de que a convergência prematura (estagnação em torno de um ótimo local) possa ser evitada ou que o problema não seja sensível à escolha desses parâmetros.

Diante do que foi apresentado, parece ser promissor o desenvolvimento de uma estratégia onde os parâmetros (probabilidade de cruzamento, taxa de perturbação e tamanho da população) do algoritmo de ED possam ser atualizados dinamicamente, ou seja, a cada geração, a partir do uso de informações sobre o processo evolutivo. Neste contexto, a próxima seção apresenta a metodologia utilizada para a atualização dinâmica dos parâmetros do algoritmo de ED.

4 ATUALIZAÇÃO DINÂMICA DOS PARÂMETROS DO ALGORITMO DE EVOLUÇÃO DIFERENCIAL

Como destacado anteriormente, uma das principais dificuldades encontradas na utilização dos métodos heurísticos é a definição dos valores dos parâmetros que serão utilizados na resolução do problema em questão. A maioria dos trabalhos apresentados na literatura ^13),(14, apesar de apresentarem resultados satisfatórios para uma grande gama de problemas, não fazem qualquer análise de sensibilidade destes parâmetros em relação à qualidade da solução obtida. Aliado a isso, as estratégias descritas por Sun e colaboradores ⁷ para a atualização do tamanho da população independem de quaisquer informações sobre o andamento do processo evolutivo, o que, do ponto de vista do processo de otimização, representa uma grande desvantagem.

Neste contexto serão apresentadas, nas próximas seções, as estratégias para a atualização dos parâmetros do algoritmo de ED. Estas serão utilizadas para o desenvolvimento do algoritmo heurístico adaptativo a que esse trabalho se propõe.

4.1 Taxa de Perturbação e da Probabilidade de Cruzamento

A atualização dinâmica da taxa de perturbação (F) e da probabilidade de cruzamento (CR) será realizada segundo a metodologia proposta por Zaharie ¹. Esta estratégia é baseada na evolução da variância da população X de tamanho NP durante o processo de otimização, conforme a equação (4.1).

(4.1)

em que

(4.2)

Segundo Zaharie ¹, o valor esperado para a variância da população, obtida após a recombinação (processo de geração de um novo candidato no algoritmo de ED), se o melhor indivíduo da população não é considerado, é dado pela equação (4.3).

(4.3)

Seja X(g) a população obtida na geração g-1 (população inicial para a geração g). Durante a geração g a população X é transformada em X' (por recombinação) e em X'' (pelo processo de seleção do melhor indivíduo).

Definindo γ como sendo Var(X(g+1))/Var(X(g)), a tendência da variância da população pode ser analisada. Se γ < 1 pode-se compensar um incremento na variância e o processo de convergência é acelerado, mas com risco de convergência prematura. Por outro lado, se γ > 1 pode-se compensar um decremento na variância, evitando assim convergência prematura. Assim sendo, a ideia é escolher os parâmetros F e CR de forma a compensar os efeitos da aplicação dos operadores de recombinação e seleção em gerações anteriores. A partir da manipulação algébrica da equação (4.3), os parâmetros F, CR e NP podem ser relacionados com a variância da população, conforme a equação (4.4).

(4.4)

em que

(4.5)

A equação (4.4) pode ser resolvida com relação ao parâmetro F, isto é:

(4.6)

com η ≡ NP(λ - 1) +2CR (2 - CR) e F _min o valor mínimo para F ∊ [0,12].

Analogamente, uma expressão para CR também pode ser obtida:

(4.7)

em que CR _min é o valor mínimo para CR ∊ [0,01 1] e o parâmetro λ é definido como apresentado na equação (4.5).

Desta forma é possível, a partir da diversidade da população, atualizar os valores de F e de CR de maneira a procurar evitar a convergência prematura, bem como buscar a região onde, provavelmente, o ótimo global está localizado.

4.2 Atualização do Tamanho da População

A atualização do tamanho da população é feita através da definição da taxa de convergência (TC) conforme a equação (4.8)

(4.8)

em que f_average é o valor médio da função objetivo da população e f_worst é o pior valor da função objetivo da população, respectivamente.

A relação acima nada mais é do que uma forma de mensurar a homogeneidade da população. No início do processo evolutivo, como os valores de f_average e de f_worst são bem distintos (diversidade da população), o valor de TC é próximo de zero. Com o aumento do número de gerações, o valor de TC tende a se tornar mais próximo da unidade, já que a população tende a se tornar homogênea.

A equação (4.9) é proposta para descrever a relação entre o tamanho da população e TC.

(4.9)

em que NP _min e NP _max representam os valores mínimo e máximo para o número de indivíduos da população, respectivamente, e o operador round indica o arredondamento para o inteiro mais próximo. Neste caso, o número de indivíduos da população pode ser atualizado dinamicamente a partir de informações sobre a convergência (ou não) do processo evolutivo.

No caso dessa estratégia, três possibilidades devem ser consideradas:

• Se o tamanho da população não se altera, nenhuma providência é tomada;

• Se o tamanho da população na próxima geração diminui com relação à geração corrente, os indivíduos que farão parte da nova população (com tamanho reduzido) são selecionados a partir da ordenação dos mesmos em relação ao valor da função objetivo, formando assim uma nova população com dimensão igual ao tamanho que foi determinado;

• Se o tamanho da população na próxima geração aumenta em relação à geração corrente, gera-se a diferença de indivíduos através do operador empregado para a geração de indivíduos no algoritmo de ED, e esses são então acrescentados à população atual paraformar a nova população.

Deve ser ressaltado que a estratégia proposta apresenta como principal vantagem o uso de informações sobre o andamento do processo evolutivo do estudo de caso em questão, para a atualização do número de indivíduos da população. Neste cenário, se a população não tende a se tornar homogênea, uma população com um número máximo de indivíduos é então considerada. Por outro lado, se a população caminha para se tornar homogênea, o algoritmo tende a utilizar uma população com um número menor de indivíduos.

Para validar a metodologia apresentada nesta seção, uma série de funções matemáticas comdiferentes graus de complexidade, bem como a obtenção de um protocolo otimizado necessário para o tratamento de tumores, serão utilizados como estudos de caso. Assim, a próxima seção apresenta a modelagem matemática referente ao fenômeno de crescimento de tumores, bem como a formulação do PCO multi-objetivo.

5 MODELAGEM MATEMÁTICA DO CRESCIMENTO DE TUMORES ATRAVÉS DE UM PCO MULTI-OBJETIVO

A modelagem matemática de sistemas biológicos é uma tarefa que não é simples devido à complexidade inerente destes sistemas. Apesar da dificuldade, vários modelos podem ser encontrados na literatura ¹⁵. Dentre estes, a formulação matemática proposta por Pillis & Radunskaya ⁹ para descrever o tratamento do tumor por quimioterapia, é uma das mais investigadas hoje em dia. Neste modelo, N, T e I representam o número de células normais, tumorais e do sistema imunológico num determinado instante de tempo t, respectivamente. N e T seguem uma lei de crescimento logístico no modelo, onde estas competem por recursos disponíveis enquanto as células do sistema imunológico e as células tumorais competem dentro do esquema presa-predador e, a variável de controle, responsável pela administração de drogas, é representada por u (concentração de drogas). Matematicamente, o modelo adimensional que descreve este fenômeno é representado pelo sistema de equações diferenciais a seguir ⁹:

(5.1)

(5.2)

(5.3)

onde os parâmetros α_i , i=1,2,3, estão relacionados com as taxas de morte celular, b_j , j=1,2, representam a capacidade de sobrevivência, c_k , k=1,2,3,4, são as taxas de competição, d ₁ representa a taxa de mortalidade, r_p , p=1,2, são as taxas de crescimento associado a cada espécie de célula, s é a taxa de fonte de entrada de células imunes, α está relacionado com a declividade inversa da curva de resposta imune e ρ é a taxa de resposta imune.

No modelo formado pelas Eqs. (5.1)-(5.3), a presença de células tumorais estimula a resposta imunológica, representada pelo termo de crescimento não-linear positivo ρIT/(α +T). Já o termo TN aproxima a probabilidade de que ocorra um encontro entre as células tumorais e normais com uma consequente disputa de recursos disponíveis e o termo IT descreve a probabilidade de encontros entre a presa (T) e o predador (I).

Neste trabalho, diferentemente do problema proposto por Pillis & Radunskaya ⁹, considera-se como objetivos a minimização do número de células tumorais e a minimização da concentração total de drogas administradas ao paciente, definidas respectivamente como:

(5.4)

(5.5)

em que t_f corresponde ao tempo total de tratamento.

Assim, o perfil otimizado da variável de controle u deve ser determinado de modo a diminuir o número de células tumorais e diminuir a concentração total de drogas administradas ao paciente ao final do tratamento.

O problema proposto é um definido como um Problema de Controle Ótimo Multi-objetivo (PCOM), cuja solução é denominada de Curva de Pareto. Esta curva representa um conjunto de soluções de compromisso para um problema multi-objetivo, diferentemente do problema com um único objetivo onde a solução ótima é caracterizada pelo melhor valor obtido. Em contraste com a otimização de um único objetivo, a solução de um problema multi-objetivo é mais um conceito do que uma definição ¹⁶. A noção de "ótimo" foi inicialmente proposta por Francis Ysidro Edgeworth ¹⁷, sendo a mesma aperfeiçoada por Vilfredo Pareto ¹⁸. Esta é baseada na convicção intuitiva de que um ponto x ^* é tomado como ótimo se "nenhum critério utilizado pode melhorar a solução, sem piorar pelo menos um outro critério". Todavia, o ótimo de Edgeworth-Pareto quase sempre não nos fornece uma solução única, mas sim um conjunto de soluçõesdenominadas não inferiores ou soluções não-dominadas ¹¹. O Postulado de Edgeworth-Pareto ou Postulado de Pareto, como é mais conhecido, tem sido base para o desenvolvimento de teoremas importantes na teoria de otimização multi-objetivos ^19),(20.

6 RESULTADOS E DISCUSSÃO

Com o intuito de se comparar o desempenho da estratégia proposta neste trabalho e denominada como Evolução Diferencial Auto-adaptativa (EDA) (com os parâmetros F, CR e NP sendo atualizados conforme descrito na Seção 4), será utilizado o algoritmo de ED considerando diferentes configurações: ED_α (CR=0,5 e F=0,3), ED_b (CR=0,5 e F=0,5), ED_c (CR=0,8 e F=1,2) e a estratégia ED_rand em que CR e F são atualizados aleatoriamente (0≤ CR≤ 1 e 0≤ F≤ 2). Para cada uma destas estratégias será considerada uma população com 50 indivíduos, estratégia 7 (ED/Rand/1/Bin) e, como critério de parada, a homogeneidade da população (módulo da diferença entre o valor médio e o pior valor da função objetivo da geração corrente menor que 1E-10). Caso esse critério de parada não seja satisfeito o algoritmo evolui até atingir 1000 gerações. Cada uma destas estratégias é executada dez vezes considerando o conjunto de sementes dado por [0 1 2 ... 9] para a inicialização do gerador de números aleatórios. O tratamento das restrições de desigualdade foi realizado através do método da penalidade estática com parâmetro de penalização da ordem de 1E+20. Para o EDA, a seguinte faixa para o tamanho da população é considerada: NP _min ≤ NP ≤ NP _max, em que NP _min = 5 e NP _max= 50. Para mensurar a redução percentual, em termos do número de avaliações da função objetivo (n_eval ) obtida com a aplicação do algoritmo EDA, será definido um parâmetro comparativo, a saber, a taxa de redução (TR):

(6.1)

em que i representa a estratégia utilizada (i=a,b,c ou rand).

A metodologia apresentada será empregada inicialmente para a resolução de funções matemáticas com diferentes graus de complexidade e, a seguir, para a resolução do PCO multi-objetivo.

6.1 Aplicações Matemáticas

Para avaliar a metodologia proposta são consideradas as seguintes funções matemáticas a serem minimizadas ^12),(21.

• Função F ₁: Mínimo global igual a -18,55, com 0 ≤ x ₁, x ₂ ≤ 10.

(6.2)

• Função F ₂: Mínimo global igual a 44, com - 100 ≤ x_i 100(i=1,...,4).

(6.3)

sujeito a:

(6.4)

(6.5)

(6.6)

• Função F ₃: Mínimo global igual a -10,8181, com 0 ≤ x_i ≤ 10(i=1,2 e 3), 0 ≤ x_i ≤ 1 (i=4 e 5) e 0 ≤ x ₆ ≤ 2.

(6.7)

sujeito a:

(6.8)

(6.9)

(6.10)

(6.11)

(6.12)

A Tabela 2 apresenta os valores médios obtidos com a aplicação das estratégias evolutivas (ED _α , ED_b , ED_c , ED_rand e EDA) nas funções matemáticas consideradas. Nesta tabela é possível observar que ambas as estratégias foram capazes de encontrar o ótimo global. Conforme esperado, a escolha dos parâmetros do algoritmo de ED influencia o número de avaliações da função objetivo (n_eval ) necessárias para a resolução dos estudos de caso propostos. Em termos do n_eval observa-se que a estratégia adaptativa proposta reduziu, em média, o número de avaliações da função objetivo em cerca de 35,77%, 40,40%, 44,66% e 41,05% em relação aos algoritmos ED _α , ED_b , ED_c e ED_rand , respectivamente. Considerando os parâmetros CR e F como sendo atualizados aleatoriamente, a probabilidade de que a população se torne homogênea aumenta. Todavia, esse procedimento não faz com que, necessariamente, se tenha a redução no n_eval para todas as estratégias evolutivas.

Tabela 2:
Resultados médios obtidos com a aplicação das estratégias evolutivas nas funções matemáticas.

F1	EDa	EDb	EDc	EDrand	EDA
x1	9,038/1E-10*	9,039/1E-11	9,039/1E-11	9,038/2E-11	9,039/2E-11
x2	8,668/2E-11	8,669/1E-11	8,668/2E-10	8,669/1E-10	8,669/1E-10
f	-18,555/1E-13	-18,554/2E-12	-18,555/4E-14	-18,555/2E-14	-18,555/2E-14
neval	4640/247	5395/227	7895/317	4565/349	997/322
TR (%)	78,51	81,51	87,37	78,15	-
F2	EDa	EDb	EDc	EDrand	EDA
x1	0,000/1E-8	0,000/1E-9	0,000/2E-8	0,000/2E-8	0,000/1E-9
x2	0,999/2E-8	0,999/1E-8	1,000/2E-9	0,999/6E-9	1,000/2E-9
x3	1,999/1E-9	2,001/2E-9	1,999/2E-8	1,999/3E-9	2,000/2E-9
x4	-1,000/3E-8	-0,999/1E-8	-1,000/1E-9	-0,999/1E-8	-0,999/2E-9
f	6,000/3E-10	6,000/1E-11	6,000/3E-9	5,999/2E-10	5,999/1E-9
neval	8750/450	9550/1080	10160/850	9250/1060	7590/590
TR (%)	13,25	20,52	25,29	17,94	-
F3	EDa	EDb	EDc	EDrand	EDA
x1	0,000/1E-9	0,000/2E-8	0,000/1E-9	0,000/5E-8	0,000/3E-8
x2	5,696/2E-9	5,697/2E-8	5,697/2E-8	5,698/1E-9	5,697/5E-8
x3	0,000/2E-9	0,000/1E-8	0,000/3E-8	0,000/1E-7	0,000/6E-9
x4	1,000/2E-9	1,000/1E-9	0,999/1E-9	0,999/5E-8	1,000/5E-8
x5	1,000/6E-8	0,999/4E-7	1,000/4E-7	0,999/1E-7	0,999/5E-7
x6	0,121/2E-8	0,122/1E-8	0,122/2E-7	0,120/2E-7	0,121/6E-8
f	-10,818/1E-9	-10,818/1E-8	-10,818/4E-8	-10,818/2E-9	-10,818/1E-8
neval	35390/12454	36980/12345	37990/9898	40980/10495	29880/10970
TR (%)	15,56	19,19	21,34	27,08	-

• ∗Desvio padrão.

A Figura 3 apresenta a taxa de convergência bem como a atualização do tamanho da população ao longo do processo evolutivo para a função F ₁. Neste caso, observa-se que a população caminha para a homogeneidade, o que implica que a taxa de convergência tende a unidade e, por consequência, o tamanho da população assume o seu valor mínimo. Assim, com aproximadamente 15 gerações, a população é praticamente homogênea, não sendo necessário trabalhar com o seu valor máximo, o que na prática implica na redução do n_eval requeridas para a resolução do problema, conforme observado na Tabela 2.

Figura 3:
Influência da taxa de convergência - TC e NP/50 (semente igual a zero).

Já a Figura 4 apresenta a atualização dos parâmetros CR e F ao longo do processo evolutivo para a função F ₁.

Figura 4:
Parâmetros do algoritmo de ED - CR e F/2 (semente igual a zero).

Ao se atualizar dinamicamente estes parâmetros a população torna-se homogênea mais rapidademente do que nas outras estratégias consideradas, conforme o n_eval apresentado na Tabela 2. Neste cenário, essa atualização faz com que a chance de existir diversidade na população aumente e, como consequência, novas regiões do espaço de projeto possam ser exploradas em busca da solução global. Este comportamento também é observado com a aplicação da estratégia ED_rand , todavia, de forma aleatória, isto é, sem que nenhuma informação sobre o processo evolutivo seja utilizada. Além disso, observa-se maior esforço computacional conforme o valor de n_eval requeridos em cada estudo de caso.

6.2 Problema de Controle Ótimo

Nesta aplicação são considerados os parâmetros propostos por Pillis & Radunskaya ^9),(10 e que não é enfocado nenhum tipo de tumor em particular, apesar de, qualitativamente o modelo representar o crescimento do tumor. Estes parâmetros foram escolhidos pelos autores de forma a garantir a estabilidade do sistema. Uma descrição detalhada sobre a escolha dos parâmetrospode ser encontrada em ^9),(10),(15.

Para um melhor entendimento da metodologia proposta, os seguintes pontos devem ser destacados:

• • O perfil da variável de controle é aproximado por elementos de controle constantes, isto é:

(6.13)

Como a droga deve ou não ser adminstrada ao paciente, isto caracteriza um típico problema on-off (se a droga é administrada ao paciente, u é igual a unidade, caso contrário é igual a zero). Assim, o PCO original pode ser reescrito como um problema de otimização que consiste na determinação dos instantes de tempo (t_i ,i = 0,...,n_elem ) em que cada ação de controle estará vigente (em que n_elem é o número de elementos de controle). Em termos práticos, a primeira ação de controle será ativada, isto é, u=1 para o intervalo de tempo (t ₀ = 0 ≤ t ≤ t ₁). Já para o segundo intervalo de tempo (t ₁ ≤ t ≤ t ₂), o controle não estará ativo (u=0). Este procedimento é repetido para todos os intervalos de tempo considerados de modo que se tenha u=[1 0 1 0 ...]. A partir de algumas execuções prévias foi verificado que o valor de n_elem a ser utilizado nesta aplicação será igual a 10. Valores acima não implicaram em mudanças significativas no valor da função objetivo;

• • Para a resolução do problema multi-objetivo foi considerado o Método da Soma Ponderada ²²:

(6.14)

em que ω (0 ≤ ω ≤ 1) é o peso associado a cada um dos objetivos. É importante ressaltar que, como o modelo é adimensionalizado e os objetivos tem a mesma ordem de grandeza (conforme simulações prévias), os valores de ω representam a importância relativa de cada um dos objetivos considerados, não sendo necessário realizar qualquer tipo de normalização/adimensionalização. Com a aplicação deste método será obtida a Curva de Pareto, que representa a um conjunto de soluções de compromisso entre os dois objetivos ^22),(11;

• • A integração do sistema de equações diferenciais é realizada a partir da aplicação do Método de Runge-Kutta-Fehlberg ²³;

• • Parâmetros adimensionais do modelo de Pillis & Radunskaya ⁹: α ₁ = 0,2; α ₂ = 0,3; α ₃ = 0,1; b ₁ = b ₂ = 1; c ₁ = c ₃ = c ₄ = 1; c ₂ = 0,5; d ₁ = 0,2; r ₁ = 1,5; r ₂ = 1 s=0,33; α =0,3; ρ =0,01 e t_f = 150. Condições iniciais: N(0)=0,9 e I(0)=T(0)=0,25 (10⁹ células).

Para fins de comparação, a Figura 5 apresenta os perfis de células sem a ação da variável de controle (u=0). Neste caso, observa-se que, se nenhuma droga é administrada ao paciente, o número de células tumorais aumenta e o número de células normais e imunes diminui, o que está em concordância com o esperado biologicamente.

Figura 5:
Perfis de células sem a ação da variável de controle.

A Tabela 3 apresenta os valores médios obtidos (para os dois objetivos) com a aplicação das estratégias evolutivas (ED_a , ED_b , ED_c , ED_rand e EDA) para a resolução do PCO.

Tabela 3:
Resultados médios obtidos com a aplicação das estratégias evolutivas no PCO.

ω	EDa		EDb		EDc		EDrand		EDA
0,10	14,01	5,19	14,01	5,19	14,01	5,19	14,00	5,19	14,00	5,19
0,20	10,97	5,70	10,97	5,70	10,75	5,75	10,97	5,70	10,97	5,70
0,30	9,10	6,30	9,22	6,25	9,10	6,30	9,00	6,34	9,10	6,30
0,40	7,80	6,99	7,78	7,01	7,76	7,02	7,78	7,01	7,76	7,02
0,50	6,79	7,81	6,79	7,81	6,72	7,88	6,79	7,81	6,77	7,83
0,60	5,91	8,89	5,91	8,89	6,03	8,76	5,91	8,89	5,91	8,89
0,70	5,22	10,18	5,22	10,18	5,22	10,18	5,22	10,18	5,22	10,18
0,80	4,59	12,11	4,59	12,11	4,59	12,11	4,59	12,11	4,59	12,11
0,90	4,02	15,48	4,02	15,46	4,01	15,52	4,02	15,48	4,02	15,48
TR (%)	46,65*		48,40		50,82		23,17		-

• ∗Valores médios.

Com os resultados apresentados na Tabela 3 conclui-se que todas as estratégias utilizadas foram capazes de encontrar um conjunto de soluções de compromisso que formam a solução do problema multi-objetivo. A estratégia adaptativa proposta reduziu em cerca de 46,65%, 48,40%, 50,82% e 23,17% o número de avaliações da função objetivo em relação aos algoritmos ED_a , ED_b , ED_c , ED_rand , respectivamente. Neste caso, houve uma redução significativa no n_eval requeridas pelo EDA e, por consequência, no tempo total de processamento, o que para aplicações mais complexas ou que exigem um tempo maior para a avaliação de cada candidato, representa um ganho relevante.

A Figura 6 apresenta a solução de compromisso obtida pelas estratégias heurísticas consideradas. Nesta é possível observar que ambas as estratégias encontraram a mesma forma para essa curva, que representa o ótimo segundo o critério de dominância proposto por Pareto. Todavia, é importante ressaltar que cada solução foi obtida com um número diferente de avaliações da função objetivo, conforme apresentado na Figura 7.

Figura 6:
Curva de Pareto.

Figura 7:
Número de avaliações da função objetivo.

Com a obtenção da Curva de Pareto considerando diferentes valores de ω, foi escolhido o ponto ω = 0,5 para a simulação dos perfis das variáveis de estado (Fig. 8) e de controle (Fig. 9). Este ponto foi escolhido por representar uma combinação onde ambos os objetivos tem a mesma importância, isto é, representa uma solução de compromisso entre a minimização do número de células tumorais e a minimização da concentração máxima de drogas administradas ao paciente. Na Figura 9 observa-se que em apenas dois intervalos de tempo a variável de controle está ativa, ou seja, a droga está sendo administrada ao paciente. Do ponto de vista clínico este resultado é bem interessante, já que a minimização da quantidade de drogas que devem ser administradas implica na minimização dos prováveis efeitos colaterias que cada paciente esta sujeito durante o tratamento. É importante ressaltar que, mesmo com essa pequena quantidade de drogas administrada (em relação a outros pontos da Curva de Pareto), o número de células tumorais diminui ao longo do tratamento, representando um protocolo otimizado para adminstração de drogas em pacientes com câncer.

Figura 8:
Perfis de células para ω = 0,5.

Figura 9:
Perfil da variável de controle para ω = 0,5.

7 CONCLUSÕES

Este trabalho objetivou o desenvolvimento de um algoritmo auto-adaptativo para a atualização dinâmica dos parâmetros do algoritmo de ED. Tal metodologia é baseada no conceito de diversidade da população e na definição de taxa de convergência para métodos heurísticos. A partir da análise dos resultados obtidos com a aplicação em funções puramente matemáticas e, também, levando em conta a resolução de um PCO multi-objetivo para fins da minimização do número de células tumorais e da minimização da quantidade de drogas que deve ser administrada ao paciente durante o tratamento, conclui-se que a utilização de informações sobre o andamento do processo evolutivo parece ser de grande importância para uma atualização adequada dos parâmetros do algoritmo de ED. Neste caso, com a abordagem proposta foi possível reduzir o número de avaliações da função objetivo requeridas para a resolução de cada estudo de caso quando comparado com o algoritmo de ED (considerando diferentes configurações para os parâmetros).

Com relação ao PCO multi-objetivo proposto, foi possível determinar um conjunto de soluções de compromisso entre os dois objetivos a partir da aplicação do Método da Soma Ponderada em associação com o algoritmo auto-adaptativo proposto. A Curva de Pareto obtida é constituída por uma série de pontos, onde cada um representa um protocolo otimizado que relaciona, de acordo com o peso considerado, os dois objetivos. Neste caso, pode-se, a partir da definição de uma estratégia própria para escolher um destes pontos, chegar a um protocolo otimizado para a administração de drogas para ser implementado na prática.

Cabe ressaltar que o objetivo deste trabalho não foi resolver um problema multi-objetivo através da aplicação de técnicas fundamentadas no critério de dominância de Pareto, mesmo sabendo já que estas são as mais apropriadas para a resolução deste tipo de problema, mas sim a apresentação de uma metodologia para a atualização dinâmica dos parâmetros do algoritmo de ED, com aplicação em um problema multi-objetivo de grande relevância nos dias atuais. Com relação a aplicação do Método da Soma Ponderada, ressalta-se que nenhuma dificuldade no que tange a normalização/adimensionalização foi observada, já que o modelo é adimensionalizado e os objetivos tem a mesma ordem de grandeza.

Finalmente, ressalta-se que os parâmetros que caracterizam o problema definido na Seção 6 não representam um tipo específico de tumor. Para essa finalidade é necessário a definição dos parâmetros que caracterizam o modelo diferencial que representa esse fenômeno.

Como sequência deste trabalho pretende-se avaliar o desempenho da metodologia proposta no projeto de sistemas robustos, ou seja, de sistemas em que a solução, dita robusta, é capaz de mitigar as incertezas a eles inerentes.

AGRADECIMENTOS

Os autores agradecem a FAPEMIG e a CAPES pelo suporte financeiro deste trabalho. O se gun do autor agradece o suporte financeiro através da FAPEMIG e do CNPq (INCT-EIE).

REFERÊNCIAS

• ¹
  D. Zaharie. Critical Values for the Control Parameters of Differential Evolution Algorithms. Proceedings of the 8^th International Conference on Soft Computing, 62-67, (2002).
• ²
  R. Storn & K. Price. Differential Evolution: A Simple and Efficient Adaptive Scheme for Global Optimization over Continuous Spaces. International Computer Science Institute, 12 (1995), 1-16.
• ³
  R. Storn, K. Price & J.A. Lampinen. Differential Evolution - A Practical Approach to Global Optimization. Springer - Natural Computing Series, (2005).
• ⁴
  F.S. Lobato, L.C. Oliveira-Lopes & V. Steffen Jr. Uma Abordagem Evolutiva Baseada em Modelos Caóticos de Busca Associado ao Algoritmo de Evolução Diferencial. 17° POSMEC - Simpósio do Programa de Pós-Graduação em Engenharia Mecânica, Uberlândia-Brasil, (2007).
• ⁵
  M.S. Tavazoei & M. Haeri. Comparison of Different One-dimensional Maps as Chaotic Search Pattern in Chaos Optimization Algorithms. Applied Mathematics and Computation, 187 (2007), 1076-1085.
• ⁶
  D. Yang, G. Li & G. Cheng. On the Efficiency of Chaos Optimization Algorithms for Global Optimization. Chaos, Solutions and Fractals, 34 (2007), 366-1375.
• ⁷
  S.Y. Sun, G. Quiang Yc, Y. Liang, Y. Liu & Q. Pan. Dynamic Population Size based Particle Swarm Optimization, ISICA 2007, Spring-Verlag Berlin, 382-392, (2007).
• ⁸
  C.C. Parker & D.P. Dearnaley. Radical Radiotherapy for Prostate Cancer. Cancer Treatment Reviews, 29 (2003), 161-169.
• ⁹
  L.G. Pillis & A. Radunskaya A. A Mathematical Tumor Model with Immune Resistance and Drug Therapy: An Optimal Control Approach. Journal of Theoretical Medicine, 3 (2001), 79-100.
• ¹⁰
  L.G. Pillis & A. Radunskaya A. The Dynamics of an Optimally Controlled Tumor Model: A Case Study. Mathematical and Computer Modelling, 37 (2003), 1221-1244.
• ¹¹
  F.S. Lobato. Otimização Multi-objetivo para o Projeto de Sistemas de Engenharia, Tese de Doutorado. Programa de Pós-Graduação em Engenharia Mecânica, Universidade Federal de Uberlândia, (2008).
• ¹²
  R.L. Haupt & S.E. Haupt. Practical Genetic Algorithms. First edition. INC. John Wiley & Sons, (1998).
• ¹³
  B.V. Babu, P.G. Chakole & J.H.S. Mubeen. Multi-objective Differential Evolution (MODE) for Optimization of Adiabatic Styrene Reactor. Chemical Engineering Science, 60 (2005), 4822-4837.
• ¹⁴
  Ji-Pyng Chiou & F.S. Wang. Hybrid Method of Evolutionary Algorithms for Static and DynamicOptimization Problems with Application to a Fed-Batch Fermentation Process. Computers and Chemical Engineering, 23 (1999), 1277-1291.
• ¹⁵
  J.C. Silveira. Modelagem Matemática do Crescimento do Câncer de Próstata com Controle Ótimo do Tratamento através da Radioterapia, Dissertação de Mestrado, (2007).
• ¹⁶
  R.T. Marler & J.S. Arora. Survey of Multi-objective Optimization Methods for Engineering. Structural Multidisciplinary Optimization, 26 (2004), 369-395.
• ¹⁷
  F.Y. Edgeworth. Mathematical Physics, London-England, First Edition, (1881).
• ¹⁸
  V. Pareto. Cours D'Economie Politique, Vol. I and II, F. Rouge, Lausanne, First Edition, (1896).
• ¹⁹
  A. Osyczka. Multicriterion Optimization in Engineering with Fortran Programs, Ellis Horwood Limited, England, First Edition, (1984).
• ²⁰
  K. Deb. Multi-Objective Optimization using Evolutionary Algorithms, John Wiley & Sons, New York, First Edition, (2001).
• ²¹
  M.M. Ali & A. Torn. Population Set Based Global Optimization Algorithms: Some Modifications and Numerical Studies. Computers and Operations Research, 31 (2004), 1703-1725.
• ²²
  G.N. Vanderplaats. Numerical Optimization Techniques for Engineering Design. Vanderplaats Research and Development, Inc., 3rd ed, (1999).
• ²³
  M.E. Davis. Numerical Methods and Modeling for Chemical Engineers. Dover Books on Mathematics, Reprint Edition, (2013).

Datas de Publicação

• Publicação nesta coleção
  May-Aug 2016

Histórico

• Recebido
  05 Jan 2016
• Aceito
  28 Jul 2016