\documentclass{TEMA}
\usepackage[brazil]{babel}      % para texto em Português
\usepackage[latin1]{inputenc}   % para acentuação em Português
\usepackage{latexsym,amssymb,amsmath,pifont}
\usepackage[latin1]{inputenc}
\usepackage[brazil]{babel}
\usepackage{memhfixc}
\usepackage{subfigure}
\usepackage[dvips]{graphicx}
\usepackage{makeidx}
\usepackage{setspace}
\usepackage{pdflscape}
\usepackage{multirow}
\usepackage{rotating}
\usepackage{longtable}
\usepackage[usenames]{color}
\usepackage{lettrine}
\usepackage[T1]{fontenc}
\usepackage{ae}
\usepackage{indentfirst}
\usepackage{yfonts}
\usepackage[all]{xy}
\usepackage{epsfig}

\newcommand{\B}{{\tt\symbol{92}}}
\newcommand{\til}{{\tt\symbol{126}}}
\newcommand{\chap}{{\tt\symbol{94}}}
\newcommand{\agud}{{\tt\symbol{13}}}
\newcommand{\crav}{{\tt\symbol{18}}}

\begin{document}

%********************************************************
\title
    {Controle Ótimo de Sistemas Algébrico-Diferenciais Chaveados usando
    o Algoritmo de Busca Gravitacional}

\author
    {A. A. PFEIFER
     \thanks{adriene$_{-}$br@yahoo.com.br}\,,
     Instituto Federal de Educação, Ciência e Tecnologia de Goiás,
     Campus Itumbiara, CEP 75524-010 - Itumbiara - GO, Brasil.
     \\\\
    %
     F. S. LOBATO
     \thanks{fslobato@feq.ufu.br. Os autores agradecem o suporte financeiro da FAPEMIG e do CNPq.}\,,
     Faculdade de Engenharia Química, Universidade Federal de Uberlândia,
     Campus Santa Mônica - CX 593, CEP 38408-100 - Uberlândia - MG, Brasil.}

\criartitulo

\runningheads {Pfeifer e Lobato}{Controle Ótimo de EADs usando o Algoritmo de Busca Gravitacional}

\begin{abstract}
{\bf Resumo}. O Problema de Controle Ótimo Chaveado (PCOC) consiste na determinação do perfil da variável de controle que minimiza uma função objetivo sujeito a restrições algébrico-diferenciais definidas por fases. Tradicionalmente, este proble- ma tem sido resolvido usando técnicas clássicas, isto é, por abordagens que fazem uso de informações sobre o gradiente da função objetivo e das restrições. Este tratamento numérico consiste na manipulação algébrica do problema original (obtenção das equações de sensibilidade) e da resolução de um problema de valor no contorno altamente dependente das estimativas iniciais. Neste contexto, o presente trabalho tem por objetivo a resolução de PCOCs usando o algoritmo de Busca Gravitacional. Tal estratégia de busca é fundamentada na lei Newtoniana de gravidade para a ge- ração de candidatos em potencial para a resolução de problemas de otimização. Os resultados obtidos são comparados com aqueles encontrados pelo Algoritmo de Levenberg-Marquardt e com uma versão híbrida.

{\bf Palavras-chave}. Controle Ótimo, Equações Algébrico-Diferenciais Chaveados, Algoritmo de Levenberg-Marquardt, Algoritmo de Busca Gravitacional, Hibridização.
\end{abstract}

%********************************************************
\newsec{Introdução}

O Problema Controle Ótimo (PCO), também conhecido como Problema de Otimização Dinâmica,
consiste na obtenção do perfil da variável de controle que minimiza um determinado índice de
desempenho (função objetivo) sujeito a restições algébrico-diferenciais. A principal dificuldade encontrada
na resolução deste proble- ma é a presença de singularidades, que resultam da
flutuação do índice diferencial ou da presença da variável de controle na forma linear na função Hamiltoniana, definida a partir
da aplicação das condições de otimalidade \cite{BrHo75,Feeh98,Xu01,Lobato,Pfeifer07}, e de restrições de desigualdade com
variáveis de estado \cite{BrHo75,Feeh98}.

Durante décadas a resolução deste problema ficou restrita ao uso
de métodos determinísticos, a saber, métodos diretos, indiretos e
híbridos \cite{LoBi89,Lobato,BrHo75,Bulirsch2,Bulirsch3,Feeh98}. Nos últimos anos, o
desenvolvimento de novos métodos de otimização, associado a
sofisticação dos computadores, possibilitou o aprimoramento de
estratégias sistemáticas para a resolução de PCO. Dentre estas estratégias, destaca-se o
Algoritmo de Busca Gravitacional (ABG) \cite{Rashedi}, devido aos bons resultados obtidos em aplicações distintas na ciência e na engenharia.
Este fundamenta-se na lei Newtoniana de gravidade para atualizar/gerar
candidatos em potencial para a solução de problemas de otimização.

No contexto da otimização dinâmica, Xu \cite{Xu01} propôs o Problema de Controle Ótimo Chaveado (PCOC) com índice de desempenho do tipo Bolza \cite{BrHo75}. Este consiste na definição de diferentes sistemas de equações algébrico-diferenciais válidas para determinadas regiões do domínio de integração. Basicamente, pode-se dizer que este problema consiste, além da determinação do perfil da variável de controle, da obtenção dos instantes onde cada um destes sistemas de equações diferenciais estão ativos. Segundo Xu \cite{Xu01}, o PCOC também pode ser obtido a partir da transcrição de um sistema chaveado geral em um PCO equivalente parametrizado pelos eventos. Esse problema equivalente tem a vantagem de solucionar o problema da presença de flutuações de índice decorrente das ativações e desativações de restrições ao longo da trajetória que provocam alterações na forma funcional das equações do problema \cite{Xu01,Xu03,XuAn04}. Na literatura, inúmeros trabalhos envolvendo a formulação matemática e a resolução do PCOC via métodos clássicos \cite{Xu01,Xu03,XuAn04,Pfeifer07} e via métodos evolutivos \cite{Long,LuusChen,Boubaker,Sakly09,Sakly11} podem ser encontrados.

Neste cenário, o presente trabalho tem por objetivo definir uma sistemática geral para a solução de PCOC usando o ABG, o Algoritmo de Levenberg-Marquardt (ALM) e um  híbrido envolvendo estas duas técnicas. Este artigo está estruturado como segue: a seção 2 apresenta a descrição de matemática do PCOC; na seção 3 é apresentada a concepção conceitual do ABG; na seção 4 são apresentadas duas aplicações. Finalmente, as conclusões e perspectivas futuras são apresentadas na seção 5.

\newsec{Formulação Matemática do PCOC}
\label{Sec2}
\setcounter{equation}{0}

Segundo Xu e Antsaklis \cite{XuAn04}, o PCOC pertence a uma classe particular de sistemas dinâmicos, nos quais os comportamentos contínuos e discretos interagem entre si, mas cujo estado contínuo não apresenta ``saltos" \;nos eventos. Neste caso, em cada domínio $t_{n}<t\leq t_{n+1}$, a dinâmica do processo segue uma diferente lei de formação, podendo implicar em descontinuidades durante a integração deste sistema. Matematicamente, o PCOC pode ser descrito como \cite{Xu01}:
%
\begin{equation}
\label{eq{pco1}:eps} \min J \triangleq \psi(\textbf{x}(t_{f}),t_{f})+\int_{t_{0}}^{t_{f}}L(\textbf{x},\textbf{u},t)dt
\end{equation}
\label{eq{pco2}:eps}
\begin{eqnarray}
{\bf{\dot x}}(t)=\left\{\begin{array}{ll} f_1(\textbf{x},\textbf{u},t)\;\;t_0\leq t \leq t_{s1} \\
f_2(\textbf{x},\textbf{u},t) \;\;t_{s1} < t \leq t_{s2} \\
\vdots \\
f_i(\textbf{x},\textbf{u},t)  \;\;t_{si-1} < t \leq t_{f}
\end{array} \right.
\end{eqnarray}
\begin{equation}
\label{eq{pco1res}:eps} s(\textbf{x},\textbf{u},t)\leq0
\end{equation}
\begin{equation}
\label{eq{pco5}:eps}\textbf{$\theta$}^{L}\leq \textbf{$\theta$}\leq \textbf{$\theta$}^{U} \;\;\; \theta=\{\textbf{x},\textbf{u},t_{si}\}\;\;\;
i=1, ..., n_{fases}
\end{equation}
onde $\textbf{x}$ é o vetor das variáveis de estado, $\textbf{u}$
é o vetor das variáveis de controle, $\psi(\textbf{x}(t_{f}),
t_{f})$ é o primeiro termo da função objetivo avaliado em
$t=t_{f}$ e $L(\textbf{x},\textbf{u},t)$ é o segundo termo do
funcional PCO, $s(\textbf{x},\textbf{u},t)$ é o vetor de
restrições algébricas e $f(\dot{\textbf{x}},\textbf{x},\textbf{u},t)$ é o vetor de
restrições diferenciais definidas pelos eventos $t_{si-1}$ e $t_{si}$ e $n_{fases}$ é o número de fases. Os sobrescritos $L$ e $U$
identificam, respectivamente, os limites inferior e superior das
variáveis.

Assim como no PCO, objetiva-se determinar o perfil do vetor de variáveis de controle e os eventos de modo a minimizar o funcional $J$. Para fins de simplificação, seja o sistema com uma única variável de estado $x$ e uma fase ($n_{fases}$=1). Introduzindo-se a nova variável de estado $T_{1}$ que corresponde ao evento $t_{s1}$ tem-se \cite{Xu01}:
%
\begin{align}
 \frac{{\partial T_{1} }}{{\partial t}} = 0 \;\;\;
 T_{1}(0) = t_{s1}
\end{align}

Introduzindo-se também uma variável temporal $\tau$, é estabelecida uma relação linear por partes entre
$t$ e $\tau$ \cite{Xu01}:
\begin{equation}
 t = \left\{ \begin{array}{l}
 t_0  + (T_{1}  - t_0 )\tau ,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\ {\rm{              }}0 \le \tau  < 1 \\
 T_{1}  + (t_f  - T_{1} )(\tau  - 1),\;\;\;\;\;{\rm{      }}1 \le \tau  \le 2 \\
 \end{array} \right.
\end{equation}

Claramente, pode-se perceber que se $\tau = 0 \Rightarrow t = t_0$, se $\tau = 1 \Rightarrow t = t_{s1}$ e se $\tau =
2 \Rightarrow t = t_f$. Introduzindo $T_{1}$ e $\tau$, o problema original pode ser reescrito pelo equivalente definido por fases:
\begin{equation}
\frac{{dx(\tau )}}{{d\tau }} = (T_{1}  - t_0 )f_1 (x,u) \;\; \tau \in [0,1)
\end{equation}
\begin{equation}
\frac{{dx(\tau )}}{{d\tau }} = (t_f  - T_{1} )f_2 (x,u) \;\; \tau \in [1,2]
\end{equation}
\begin{equation}
\frac{{dT_{1} }}{{d\tau }} = 0
\end{equation}

Deseja-se encontrar $t_{s1}$ e $u(\tau)$ de forma a minimizar o funcional:
\begin{gather}
\label{eq{3c}:eps}
J =\psi (x(2))+\int\limits_0^1 {(T_{1}  - t_0 )L(x,u)d\tau }+\int\limits_1^2 {(t_f-T_{1} )L(x,u)d
\tau}
\end{gather}

Definindo as funções $f$, $L$ e \textit{Hamiltoniano} ($H(x,p,u,t_{s1})$) na forma parametrizada:
\begin{gather}
\tilde f_1 (x,u,t_{s1} ) \buildrel \Delta \over = (T_{1}  - t_0 )f_1 (x,u)\;\;\;
\tilde f_2 (x,u,t_{s1} ) \buildrel \Delta \over = (t_f  - T_{1} )f_2 (x,u)\\
\tilde L_1 (x,u,t_{s1} ) \buildrel \Delta \over = (T_{1}  - t_0 )L(x,u)\;\;\;
\tilde L_2 (x,u,t_{s1} ) \buildrel \Delta \over = (t_f  - T_{1} )L(x,u)\\
H(x,p,u,t_{s1}) \buildrel \Delta \over = \left\{ \begin{array}{l}
 \tilde L_1 (x,u,T_{1} ) + p^T \tilde f_1 (x,u,T_{1} ),\ \ \tau  \in [0,1) \\
 \tilde L_2 (x,u,T_{1} ) + p^T \tilde f_2 (x,u,T_{1} ),\ \ \tau  \in [1,2] \\
 \end{array} \right.
\end{gather}
onde $p$ é a variável de co-estado (adjunta).

Aplicando as condições necessárias \cite{BrHo75}, as equações de estado, co-estado e de estacionariedade e as condições de contorno e continuidade, são obtidas:
\begin{eqnarray}
\label{eq{4c1}:eps} \frac{{\partial x}}{{\partial \tau }} = \left( {\frac{{\partial H}}{{\partial p}}}
\right)^T  &=& \tilde f_k (x,u,T_{1} )\\
\label{eq{4c2}:eps} \frac{{\partial p}}{{\partial \tau }} =  - \left( {\frac{{\partial H}}{{\partial
x}}} \right)^T  &=&  - \left( {\frac{{\partial \tilde f_k }}{{\partial x}}} \right)^T p - \left(
{\frac{{\partial \tilde L_k }}{{\partial x}}} \right)^T \\
\label{eq{4c3}:eps} 0 = \left( {\frac{{\partial H}}{{\partial u}}} \right)^T  &=& \left( {\frac
{{\partial \tilde f_k }}{{\partial u}}} \right)^T p - \left( {\frac{{\partial \tilde L_k }}{{\partial
u}}} \right)^T \\
\label{eq{4c4}:eps} x(0,T_{1} ) &=& x_0 \\
\label{eq{4c5}:eps} p(2,T_{1} ) &=& \left( {\frac{{\partial \Psi }}{{\partial x}}(x(2,T_{1} ))}
\right)^T\\
\label{eq{4c6}:eps} p(1^{-},T_{1} ) &=& p(1^{+},T_{1} )
\end{eqnarray}

As equações acima formam um sistema de EADs de valor no contorno
parame- trizado pelo evento $t_{s1}$. Deste modo, o valor ótimo do ``novo"\; funcional ($J_1$), isto é, como uma função dos parâmetros $T_{1}$ é
dado pela seguinte equação:
\begin{gather}
\label{eq{d1}:eps}
J_1 (T_{1} ) = \psi (x(2,T_{1} )) + \int\limits_0^1 {\tilde L_1 (x,u,T_{1} )d\tau }  + \int
\limits_1^2 {\tilde L_2 (x,u,T_{1} )d\tau }
\end{gather}

Diferenciando $J_1$ em relação a $T_{1}$, faz-se necessário obter as funções $\frac{\partial }{{\partial \tau }}\left( {\frac{{\partial x}}
{{\partial T_{1} }}}\right) $ e $\frac{\partial }{{\partial \tau }}\left( {\frac{{\partial u}}
{{\partial T_{1} }}}\right) $ para encontrar o valor $\frac{dJ_1}{dT_{1}}$. Assim, diferencia-se as equações
de estado, co-estado, estacionariedade e as condições de contorno e de continuidade com em relação a $T_{1}$, gerando um sistema para cada intervalo de $\tau$. Neste caso, a partir de uma estimativa inicial para $t_{s1}$ pode-se resolver o problema de otimização formado pelas equações de otimalidade associado as equações de sensibilidade. A resolução deste PCOC representa um grande desafio por demandar o conhecimento da sequência e do número de ativações e desativações ao longo da trajetória e de uma boa estimativa inicial \cite{Xu01}. Para facilitar a geração do sistema de equações de valor no contorno, Pfeifer \cite{Pfeifer07} desenvolveu uma ferramenta simbólica implementado em ambiente Maple$^{\circledR}$.

\newsec{O Algoritmo de Busca Gravitacional}
\label{Sec3}
\setcounter{equation}{0}

%\subsection{Inspiração Conceitual}
%
%A gravidade pode ser entendida como a ``tendência" que as massas têm para acelerar-se em direção umas das outras, isto é, cada partícula no universo atrai e é atraída por outra partícula \cite{Schutz}.  A forma como a força gravitacional se comporta denomina-se ``ação a distância". Isto significa que a gravidade age entre as partículas separadas, sem qualquer intermediário e sem qualquer atraso \cite{Holliday}. Fisicamente, sabe-se que a força gravitacional entre duas partículas é diretamente proporcional ao produto das suas massas e inversamente proporcional ao quadrado da distância entre elas. Já a segunda lei de Newton diz que quando uma força é aplicada a uma partícula, a sua aceleração depende da força e sua massa.
%
%Baseando-se nestas duas leis da física observa-se que existe uma força (gravidade) de interação atuando em todas as partículas do universo como ilustrado na Figura \ref{f:gravita}. Nesta figura, $F_{1j}$ é a força que atua em $M_1$ de $M_j$ e $F_1$ é a força total que atua sobre $M_1$ e provoca a aceleração do vetor $a_1$.
%
%%\begin{figure}[h]
%%\centering
%%\centerline{\includegraphics[angle=0,scale=0.4]{figuras/gravita.pdx}}
%%\caption{Resultante da força $F_1$ atuando nas massas $j$ \cite{Rashedi}.}\label{f:gravita}
%%\end{figure}
%
%Considerando os aspectos acima mencionados, podemos reescrever as leis de Newton para a nossa finalidade. A força gravitacional, $f_{ij}$, que age sobre a massa $i$ por ação da massa $j$, é proporcional ao produto das massas gravitacionais $j$ e $i$, e inversamente proporcional ao quadrado da distância entre eles. $a_i$ é proporcional a $f_{ij}$ e inversamente proporcional à massa de inércia $i$. Neste caso, pode-se reescrever as leis da física citadas anteriormente como segue:
%
%\begin{equation}
%\label{eq{desc3}:eps} F_{ij}=G\frac{M_{aj}M_{pi}}{R^2}
%\end{equation}
%\begin{equation}
%\label{eq{desc4}:eps} a_{i}=\frac{F_{ij}}{M_{ii}}
%\end{equation}
%onde $M_{aj}$ e $M_{pi}$ representam as massas gravitacionais das partículas $i$ e $j$, respectivamente, e $M_{ii}$ representa a inércia da partícula $i$.
%
%A próxima seção apresenta a descrição matemática do ABG \cite{Rashedi}.
%
%\subsection{O Algoritmo}

O ABG pertence a classe de algoritmos não-determinísticos baseados em uma população de candidatos para a resolução do problema de otimização. Este algoritmo, assim como acontece em outras abordagens evolutivas, caracteriza-se por gerar uma população inicial (aleatória) e definir uma estratégia para atualizar a população de candidatos.

No ABG o desenvolvimento dessa estratégia é fundamentada nas leis de gravidade e de movimento \cite{Rashedi}. Neste contexto, os ``agentes" \; são considerados como obje- tos e seu desempenho é mensurado por suas massas. Toda massa atrai e é atraída por outras massas  devido a força gravitacional,
causando o movimento global de todo o conjunto (população) considerado. Assim, as massas cooperaram usando uma forma direta de ``comunicação", a força gravitacional. As massas ``pesadas"\, - o que corresponde a boas soluções - movem-se mais lentamente do que as mais leves, o que garante a etapa de exploração do algoritmo.

No ABG, cada massa tem quatro especificações \cite{Rashedi}: posição, massa inercial, massa gravitacional ativa, e massa gravitacional
passiva. A posição da massa corresponde a uma solução candidata do problema de otimização, e suas massas gravi- tacionais e inerciais são
determinadas usando a função objetivo.

Seja o sistema com agentes $N$ (massas). A posição do $i$-ésimo agente é dada:
%
\begin{equation}
\label{eq{algo1}:eps} x_i=(x_{i}^{1}, ..., x_{i}^{d}, ..., x_{i}^{n}),\;\;i=1,...,N
\end{equation}
onde $x_{i}^{d}$ é a posição do $i$-ésimo agente na $d$-ésima dimensão.

Na geração $k$, define-se a força atuando na massa $i$ devido a ação da massa $j$:
%
\begin{equation}
\label{eq{algo2}:eps} F_{ij}^{d}(k)=G(k)\frac{M_{pi}(k)\times M_{aj}(k)}{R_{ij}(k)}(x_{j}^{d}(k)-x_{i}^{d}(k))
\end{equation}
onde $M_{aj}$ é a massa gravitacional ativa relativa ao agente $j$, $M_{pi}$ é a massa gravi- tacional passiva relativa ao agente $i$, $G(k)$ é a constante gravitacional na iteração $k$, e $R_{ij}(k)$ é a distância Euclidiana entre os agentes $i$ e $j$.

Para introduzir a característica estocástica ao algoritmo, Rashedi \cite{Rashedi} supôs que a força total que atua no agente $i$ é aleatoriamente ponderada pelos $d$-ésimos componentes das forças exercidas por outros agentes:
%
\begin{equation}
\label{eq{algo3}:eps} F_{i}^{d}(k)=\sum_{j=1,j\neq i}^{N}rand_j F_{ij}^{d}(k)
\end{equation}
onde $rand_j$ é um número aleatório dentro do intervalo [0,\;1].

Deste modo, segundo a lei do movimento, a aceleração do agente $i$ na iteração $k$, e na $d$-ésima direção, $a_{i}^{d}(k)$, é calculado pela seguinte equação:
%
\begin{equation}
\label{eq{algo4}:eps} a_{i}^{d}(k)=\frac{F_{i}^{d}(k)}{M_{ii}(k)}
\end{equation}
onde $M_{ii}$ é a massa inercial do $i$-ésimo agente.

Além disso, a atualização da velocidade de um agente é considerada como sendo uma fração de sua velocidade atual adicionado a sua aceleração. Portanto, a sua posição pode ser atualizada da seguinte forma:
%
\begin{equation}
\label{eq{algo6}:eps} x_{i}^{d}(k+1)=x_{i}^{d}(k)+rand_i\times v_{i}^{d}(k)+a_{i}^{d}(k)
\end{equation}
onde $rand_i$ é um número aleatório dentro do intervalo [0,\;1].

A constante gravitacional ($G(k)$) é reduzida ao longo do processo de evolutivo de acordo com a Eq. \ref{eq{algo7}:eps}:
%
\begin{equation}
\label{eq{algo7}:eps} G(k)=G_{\circ}\exp(-\alpha k/max_{iter})
\end{equation}
onde $G_{\circ}$ é a constante gravitacional inicial, $\alpha$ é um parâmetro de ajuste, $k$ é a iteração corrente e $max_{iter}$ é número máximo de iterações.

As massas gravitacional e inercial são calculadas através da avaliação da função objetivo. Uma massa mais ``pesada"\, significa um agente mais eficiente. Isto implica que os melhores agentes são os mais atraídos. Assumindo a igualdade das massas gravitacional e inercial ($M_{ai}=M_{pi}=M_{ii}=M_{i},\;\;i=1, 2, ..., N$), os valores das mesmas são calculadas usando as seguintes relações:
%
\begin{equation}
\label{eq{algo9}:eps} m_{i}(k)=\frac{f_{i}(k)-worst(k)}{best(k)-worst(k)}
\end{equation}
%
\begin{equation}
\label{eq{algo10}:eps} M_{i}(k)=\frac{m_{i}(k)}{\sum_{j=1}^{N}m_{j}(k)}
\end{equation}
onde $f_i(k)$ é a avaliação da função objetivo relativo ao agente $i$ na iteração $k$, $worst(k)$ e $best(k)$ são definidos como sendo o pior e o melhor valor da função objetivo na iteração corrente, respectivamente.

O cálculo da força $F_{i}^{d}(k)$ é realizado segundo a expressão a seguir:
%
\begin{equation}
\label{eq{algo11}:eps} F_{i}^{d}(k)=\sum_{j\in K_{best},j\neq i}rand_j\times F_{ij}^{d}(k)
\end{equation}
onde $K_{best}$ é o conjunto dos primeiros $K$ agentes com os melhores valores em termos da função objetivo. Segundo Rashedi \cite{Rashedi}, esta equação representa uma forma de se estabelecer um bom compromisso entre exploração e refinamento, de forma a reduzir a chance de ser obter um ótimo local.

É importante ressaltar que a metodologia apresentada foi proposta por Rashedi \cite{Rashedi}. Neste caso, outras formas para atualizar o valor da constante gravitacional e de modificações das massas podem ser encontradas em \cite{Nobahari11,Chatterjee}.

Na literatura pode-se encontrar algumas aplicações do ABG, dentre as quais pode-se citar:
otimização de funções matemáticas com diferentes graus de complexidade \cite{Rashedi}, desenvolvimento de um algoritmo multi-objetivo fundamentado no
operador de ordenamento por \textit{ranking} e no algoritmo ABG \cite{Nobahari11}, estudo compara- tivo envolvendo o ABG e um algorimo de enxame de partículas modificado \cite{Chatterjee}, aplicação em problemas de potência \cite{Duman} e um híbrido entre o ABG e redes neurais para resolver a equação de Wessinger \cite{Ghalambaz}.

\newsec{Resultados e Discussão}
\label{Sec4}
\setcounter{equation}{0}

Nesta seção são apresentados dois estudos de caso para a validação da metodologia proposta. Para essa finalidade alguns pontos devem ser destacados:

\begin{itemize}
  \item Algoritmo de Levenberg-Marquardt (ALM) com estimativa inicial apresentada nas tabelas a seguir. Deve ser ressaltado que o problema de valor no contorno foi resolvido usando o Método de Colocação Normal com 10 pontos e tolerância absoluta de 10$^{-6}$.
  \item Parâmetros usados no ABG \cite{Rashedi}: população com 20 agentes (gerada aleatoriamente dentro do espaço de projeto), número máximo de gerações igual a 2000, $G_{\circ}$ e $\alpha$ iguais a 100 e 20, respectivamente. Todos os estudos de caso foram executados 20 vezes em um microcomputador Intel(R) Core i5 com 3 Gb de memória para a obtenção dos resultados. Foi utilizada aproximação linear para a variável de controle definida por elementos.
  \item Método Híbrido: inicialmente executa-se o ABG (50 gerações) com posterior refinamento pelo ALM.
\end{itemize}

\subsection{Estudo de Caso 1}

Seja o PCOC descrito pelos seguintes sub-sistemas ($ss_i$) \cite{Xu01,Pfeifer07}:
%
\begin{equation}
\label{eq{resul1}} \min J=\frac{1}{2}(x_1(t_f)-e^2)^2+\frac{1}{2}(x_2(t_f)-e^2)^2+\frac{1}{2}\int_{0}^{t_f}u^2dt
\end{equation}
\begin{eqnarray}
ss_1 (0\leq t \leq ts_1)=\left\{\begin{array}{ll}  \dot{x}_1=-x_1+2x_1 \\
\dot{x}_2=x_2+u x_2
\end{array} \right. \nonumber\\
%
\label {eq{resul2}} ss_2 (ts_1 < t \leq ts_2)=\left\{\begin{array}{ll}  \dot{x}_1=x_1-3x_1 u \\
\dot{x}_2=2x_2-2u x_2
\end{array} \right.
\\
%
ss_3 (ts_2 < t \leq t_f)=\left\{\begin{array}{ll}  \dot{x}_1=2x_1+x_1 u \\
\dot{x}_2=-x_2+3u x_2
\end{array} \right.\nonumber
\end{eqnarray}
onde $x(0)$=[1 1]$^{T}$, e os eventos, $t_{si}$ ($i$=1,2), são definidos como $0 \leq t_{s1} \leq t_{s2} \leq t_f=3$. Para a resolução deste
problema, o seguinte espaço de projeto é considerado: 0 $\leq t_{si} \leq$ 2 ($i$=1,2), cuja solução global é \cite{Xu01}: $t_{s1}$=1, $t_{s2}$=2, $J$=0 e $u$=0.

A Tabela \ref{tabela01} apresenta os resultados obtidos para os eventos e a função objetivo via ABG, ALM e Híbrida (ABG-ALM). Nesta tabela, é possível observar que, apesar do menor número de iterações ($N_{iter}$) requeridas pelo ALM, existe uma grande dependência com relação a escolha da estimativa inicial. Em termos do ABG observa-se que boas estimativas são encontradas quando avalia-se o valor médio e o desvio padrão dos resultados. Já na abordagem Híbrida, como esperado, o número de gerações reduz-se consideravelmente em comparação com o ABG. Isto justifica-se pelo acoplamento de duas estratégias, uma de busca global seguida do uso de uma outra de busca local (refinamento da solução). É importante ressaltar que Pfeifer \cite{Pfeifer07} resolveu este mesmo estudo de caso via método indireto, i.e., através da geração das condições de otimalidade, encontrando função objetivo nula. Neste caso, observa-se que a abordagem indireta foi capaz obter um valor mais preciso do que os encontrados neste trabalho, justificado por sua maior precisão \cite{BrHo75}. Em termos do tempo computacional médio, observa-se que o melhor desempenho foi obtido pela abordagem Clássica (12 segundos), seguida da Híbrida (171 segundos) e do ABG (862 segundos), valores diretamente proporcionais ao número de avaliações da função objetivo requeridas por cada abordagem.
%
\begin{table}[h!]
\caption{Resultados do Estudo de Caso 1.} \label{tabela01}
\vspace{-0.65cm}
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
%
Método & $t_{s2}^{\ast}$/$t_{s2}^{\ast}$ & $t_{s1}$ & $t_{s2}$ & Solução$^{\diamond}$ & $J$ (Eq. (\ref{eq{resul1}})) & $N_{iter}$\\
\hline \hline
Xu \cite{Xu01} & 1,1/2,1 & 0,9914 & 2,0140 & C & 2,692$\times$10$^{-4}$ & 20\\\cline{1-7}
         & 0,01/1,5 & -   & -  & NC & - & - \\\cline{2-7}
         & 0,05/2,5 & -   & -  & NC & - & - \\\cline{2-7}
ALM & 0,5/1,5 &  0,99999  &  2,00000 & C & 1,155$\times$10$^{-9}$ & 33 \\\cline{2-7}
         & 0,2/1,2 & -   & -  & NC & - & - \\\cline{2-7}
         & 1,1/2,1 &  0,99999  &  2,00000 & C & 1,967$\times$10$^{-9}$ & 28 \\
         \hline
\multirow{2}{*}{ABG} & \multirow{2}{*}{\textit{rand}$^{\ast\ast}$} & 1,00068 & 2,00002 & \multirow{2}{*}{C} & 8,441$\times$10$^{-7}$ & 2000\\
         & & (1,45$\times$10$^{-7}$)$^{\dagger}$ & (2,09$\times$10$^{-7}$)$^{\dagger}$ & & (1,89$\times$10$^{-8}$)$^{\dagger}$ & - \\\cline{1-7}
         & 1,06/2,03 & 0,99999  & 2,00000  & C & 1,129$\times$10$^{-9}$ & 200+12 \\\cline{2-7}
Híbrido  & 1,09/1,93 & 0,99999 & 2,00000 & C & 2,099$\times$10$^{-9}$ & 400+10 \\\cline{2-7}
         & 0,96/2,05 & 0,99999 & 2,00000 & C & 1,934$\times$10$^{-10}$ & 600+09 \\\hline
\multicolumn{7}{l}{\footnotesize $^{\ast}$ Estimativas iniciais utilizadas no ALM e no ABG-ALM. \footnotesize $^{\dagger}$ Desvio padrão.}\\
\multicolumn{7}{l}{\footnotesize $^{\ast\ast}$ Estimativa inicial gerada aleatoriamente. \footnotesize $^{\diamond}$ C: Convergiu, NC: Não Convergiu.}\\
\end{tabular}
\end{center}
\end{table}
%

A Figura \ref{fig1} apresenta os perfis das variáveis de estado e de controle para o estudo de caso 1 usando o ALM e o ABG.
\begin{figure}[!htb]
\centerline{ \subfigure[plano de fases $x_1 \times x_2$.]
{\includegraphics[scale=0.2]{figuras/estado1.eps} } \hfil
\subfigure[controle ($u$).] {\includegraphics[scale=0.2]{figuras/controle1.eps}}}
\caption{Perfis das variáveis de estado e de controle para o estudo de caso 1.}
\label{fig1}
\end{figure}

Nestas figuras percebe-se boa concordância entre os resultados obtidos pelo ABG e pelo ALM. No caso do perfil de controle, o ABG encontrou perfil de controle ótimo constante e com valor nulo, o que não foi constatado via ALM, onde três fases distintas foram obtidas. Apesar desta diferença, oriunda da precisão adotada por cada metodologia, ambos os resultados estão em concordância com os obtidos por Xu \cite{Xu01}. Na Figura \ref{focaso1} é apresentada a evolução em termos do valor da função objetivo (normalizada) versus número de iterações para cada um dos métodos utilizados.

\begin{figure}[h]
\centering
\centerline{\includegraphics[angle=0,scale=0.25]{figuras/focaso1.eps}}
\caption{Função objetivo versus número de iterações para o estudo de caso 1 (ALM com estimativas iniciais $t_{s1}$=1,1 e $t_{s2}$=2,1; ABG-ALM com estimativas iniciais $t_{s1}$=1,09 e $t_{s2}$=1,93).}\label{focaso1}
\end{figure}
%

\subsection{Estudo de Caso 2}

O segundo estudo de caso é composto pela seguinte função objetivo, restrita por três sub-sistemas \cite{Xu01}:
%
\begin{equation}
\label{eq{resula1}} \min J=(x_1(t_f)+4,1337)^2+(x_2(t_f)-9,3569)^2+\frac{1}{2}\int_{0}^{t_f}u^2dt
\end{equation}
\begin{eqnarray}
ss_1 (0\leq t \leq ts_1)=\left\{\begin{array}{ll}  \dot{x}_1=-2x_1+u \\
\dot{x}_2=-x_2
\end{array} \right. \nonumber\\
%
\label {eq{resula2}} ss_2 (ts_1 < t \leq ts_2)=\left\{\begin{array}{ll}
\dot{x}_1=0,5x_1+5,3x_2+u \\
\dot{x}_2=-5,3x_1+0,5x_2-u
\end{array} \right.
\\
%
ss_3 (ts_2 < t \leq t_f=3)=\left\{\begin{array}{ll}
\dot{x}_1=x_1 \\
\dot{x}_2=1,5x_2+u
\end{array} \right.\nonumber
\end{eqnarray}
onde $x(0)$=[4 4]$^{T}$. Para a resolução deste problema considera-se o seguinte espaço de projeto: 0 $\leq t_{si} \leq$ 2 ($i$=1,2), cujo ótimo global é \cite{Xu01}: $t_{s1}$=1, $t_{s2}$=2, $J$=0 e $u$=0. A Tabela \ref{tabela02} apresenta os resultados obtidos para o estudo de caso 2. Assim como no caso anterior, pode-se observar que a ABG sempre convergiu para a solução, diferentemente do que acontece com o ALM. No caso da abordagem Híbrida, o número total de gerações sempre é menor do que requerido pelo ABG, resultando em uma estratégia por assim dizer robusta. Já em termos do tempo computacional médio, observa-se, assim como no caso anterior, que o melhor desempenho foi obtido pela abordagem Clássica (24 segundos), seguida da Híbrida (192 segundos) e do ABG (902 segundos), proporcionais ao número de avaliações da função objetivo requeridas por cada abordagem.

\begin{table}[h!]
\caption{Resultados do Estudo de Caso 2.} \label{tabela02}
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
%
Método & $t_{s2}^{\ast}$/$t_{s2}^{\ast}$ & $t_{s1}$ & $t_{s2}$ & Solução$^{\diamond}$ & $J$ (Eq. (\ref{eq{resula1}})) & $N_{iter}$\\
\hline \hline
Xu \cite{Xu01} & 0,8/1,8 & 1,0002 & 2,0008 & C & 6,3146$\times$10$^{-5}$ & 43\\\cline{1-7}
               & 0,8/1,8 & 1,00000 & 2,00000 & C & 2,6876$\times$10$^{-7}$ & 52\\\cline{2-7}
               & 0,9/2,1 & 1,00000 & 2,00000 & C & 3,6323$\times$10$^{-7}$ & 55\\\cline{2-7}
ALM       & 0,5/2,5 & -       & -       & NC& - & -\\\cline{2-7}
               & 0,6/2,2 & -       & -       & NC& - & - \\\cline{2-7}
               & 1,2/2,5 & -       & -       & NC& - & - \\
         \hline
\multirow{2}{*}{ABG}      & \multirow{2}{*}{\textit{rand}$^{\ast\ast}$} & 1,00002 & 2,00003 & C & 1,3334$\times$10$^{-7}$ & 2000 \\
         & & (1,55$\times$10$^{-7}$)$^{\dagger}$ & (4,99$\times$10$^{-7}$)$^{\dagger}$ & & (2,09$\times$10$^{-7}$)$^{\dagger}$ & - \\\cline{1-7}
         & 1,08/2,09 & 1,00000  & 2,00000  & C & 2,998$\times$10$^{-9}$ & 200+14 \\\cline{2-7}
Híbrido  & 1,10/1,94 & 0,99999 & 2,00000 & C & 1,198$\times$10$^{-9}$ & 400+13 \\\cline{2-7}
         & 0,97/1,99 & 1,00000 & 2,00000 & C & 1,487$\times$10$^{-10}$ & 600+14 \\\hline
\multicolumn{7}{l}{\footnotesize $^{\ast}$ Estimativas iniciais utilizadas no ALM e no ABG-ALM. \footnotesize $^{\dagger}$ Desvio padrão.}\\
\multicolumn{7}{l}{\footnotesize $^{\ast\ast}$ Estimativa inicial gerada aleatoriamente. \footnotesize $^{\diamond}$ C: Convergiu, NC: Não Convergiu.}\\
\end{tabular}
\end{center}
\end{table}
%
A Figura \ref{fig2} apresenta os perfis das variáveis de estado e de controle para o estudo de caso 2.
\begin{figure}[!htb]
\centerline{ \subfigure[plano de fases $x_1 \times x_2$.]
{\includegraphics[scale=0.2]{figuras/estado2.eps} } \hfil
\subfigure[controle ($u$).] {\includegraphics[scale=0.2]{figuras/controle2.eps}}}
\caption{Perfis das variáveis de estado e de controle para o estudo de caso 2.}
\label{fig2}
\end{figure}

Assim como observado para o estudo de caso anterior, percebe-se que apesar da sensível diferença entre os perfis de controle calculados pelo ABG e pelo ALM, relacionados com as precisões adotadas em cada metodologia, ambos os resultados também estão em concordância com os obtidos por Xu \cite{Xu01}. É importante ressaltar que o perfil de controle apresentado neste trabalho se difere do obtido por Sakly e colaboradores \cite{Sakly11}. Neste caso, estes autores obtiveram o perfil da variável de controle da ordem de 10$^{-3}$, diferentemente do ABG e do ALM que encontraram da ordem de 10$^{-5}$, mais próximo ao valor analítico ($u$=0). Na Figura \ref{focaso2} é apresentada a evolução em termos do valor da função objetivo (normalizada) versus número de iterações para cada um dos métodos utilizados.

\begin{figure}[h]
\centering
\centerline{\includegraphics[angle=0,scale=0.22]{figuras/focaso2.eps}}
\caption{Função objetivo versus número de iterações para o estudo de caso 2 (ALM com estimativas iniciais $t_{s1}$=0,9 e $t_{s2}$=2,1; ABG-ALM com estimativas iniciais $t_{s1}$=1,10 e $t_{s2}$=1,94).}\label{focaso2}
\end{figure}
%

\newsec{Conclusões}
\label{Sec5}

O presente trabalho propôs uma sistemática para a resolução de PCOCs usando o ABG e sua hibridização com uma técnica clássica fundamentada na geração de equações de sensibilidade, a saber, o ALM. Os resultados obtidos com o ABG demonstram que a abordagem empregada configura-se como uma interessante alternativa quando comparada com o ALM, apesar de demandar mais avaliações da função objetivo e, consequentemente, de um maior tempo de processamento. Do ponto de vista da hibridização, percebe-se que, devido a sensibilidade observada com a escolha da estimativa inicial (veja as Tabelas \ref{tabela01} e \ref{tabela02}), a utilização do ABG como estratégia inicial de busca, associada ao ALM, promoveu, para os estudos de caso analisados, uma estratégia que se mostrou robusta para a resolução de PCOCs. Como continuidade deste trabalho pretende-se aplicar a metodologia proposta em estudos de caso realísticos, com ênfase em problemas biotecnológicos.

\begin{abstract}
{\bf Abstract} Optimal Control Problem of Switched Systems (OCPSS) consists in determining the control variable profile that minimizes a given objective function subject to differential-algebraic constraints defined by phases. Traditionally, this problem has been solved using classical techniques through algebraic manipulation of  original problem (obtaining the sensitivity equations) and resolution, at each step, a boundary value problem highly dependent of initial estimates. In this context, the present work aim to resolve the OCPSS using the Gravitational Search Algorithm. This search strategy is based on the Newtonian law of gravity for the generation of potential candidates to solve optimization problems. The results are compared with those obtained by the Levenberg-Marquardt Algorithm and with a hybrid version.
\end{abstract}

\begin{thebibliography}{99}

\bibitem{Boubaker}
S. Boubaker and F. M'sahli, Solutions Based on Particle Swarm Optimization for Optimal Control Problems of Hybrid Autonomous Switched Systems.
{\em International Journal of Intelligent Control and Systems}, {\bf 13}(2), 128-135, (2008).

\bibitem{BrCP96}
K. E. Brenan, S. L. Campbell and L. R. Petzold, ``Numerical Solution of Initial Value Problems in Differential Algebraic Equations'',
Classics in Applied Mathematics, SIAM Philadelphia, 1996.
%
\bibitem{BrHo75}
A. E. Bryson and Y. C. Ho, ``Applied Optimal Control'', Hemisphere Publishing,Washington, 1975.
%
\bibitem{Bulirsch2}
R. Bulirsch, F. Montrone and H. J. Pesch, Abort Landing in the Presence of a Windshear as a minimax Optimal Control Problem, Part II: Multiple Shooting and Homotopy. {\em Journal of Optimization Theory and Applications}, {\bf 70}, 223-254, (1991b).
%
\bibitem{Bulirsch3}
R. Bulirsch, E. Nerz, H. J. Pesch, Combining Direct and Indirect Methods in Optimal Control: Range Maximization of a Hang Glider. {\em International Series of Numerical Mathematics}, {\bf 111}, 273-288, (1993).
%
\bibitem{Chatterjee}
A. Chatterjee and G. K. Mahanti, Comparative Performance of Gravitational Search Algorithm and
Modified Particle Swarm Optimization Algorithm for Synthesis of Thinned Scanned Concentric Ring Array Antenna,
{\em Progress in Eletromagnetics Research}, {\bf 25}, 331-348, (2010).
%
\bibitem{Duman}
S. Duman, U. Güvenç and N. Yörükeren, Gravitational Search Algorithm for Economic Dispatch
with Valve-Point Effects, {\em International Review of Electrical Engineering}, {\bf 5}, (6), 2010.
%
\bibitem{EgWD03}
M. Egerstedt, Y. Wardi and F. Delmotte, Optimal Control of Switching Times in Switched Dynamical Systems,
{\em 42nd IEEE Conference on Decision and Control} (2003).
%
\bibitem{Feeh98}
W. F. Feehery, ``Dynamic Optimization with Path Constraints'', Thesis, Massachusetts Institute of Technology, 1998.
%
\bibitem{Ghalambaz}
M. Ghalambaz, A. R. Noghrehabadi, M. A. Behrang, E. Assareh, A. Ghanbarzadeh and N.Hedayat,
A Hybrid Neural Network and Gravitational Search Algorithm (HNNGSA) Method to Solve well known Wessinger's Equation,
{\em World Academy of Science, Engineering and Technology}, {\bf 73}, (2011).
%
\bibitem{Holliday}
D. Holliday, R. Resnick, J. Walker, {\em Fundamentals of physics}, John Wiley and Sons, (1993).
%
\bibitem{Lobato}
F. S. Lobato, ``Abordagem Mista para Problemas de Otimização Dinâmica'', Dissertação de Mestrado, Faculdade de Engenharia Química,
Universidade Federal de Uberlândia, 2004.
%
\bibitem{LoBi89}
J. S. Logsdon and L. T. Biegler, Accurate Solution of Diferential-Algebraic Optimization Problems, {\em Ind. Eng. Chem. Res.}, {\bf 28}, 89-101,
(1989).
%
\bibitem{Long}
R. Long, J. Fu and L. Zhang, Optimal Control of Switched System Based on Neural Network Optimization, {\em ICIC 08 Proceedings of the 4th International
Conference on Intelligent Computing: Advanced Intelligent Computing Theories and Applications}, 799-806, (2008).
%
\bibitem{LuusChen}
R. Luus and Y. Q. Chen, Optimal Switching Control via Direct Search Optimization, {\em Proceedings of the 2003 IEEE International Symposium on Intelligent Control}, Houston, Texas, 302-306, (2003).
%
\bibitem{Nobahari11}
H. Nobahari, M. Nikusokhan and P. Siarry, Non-dominated Sorting Gravitational Search Algorithm,
{\em ICSI 2011: International conference on swarm intelligence}, Cergy, France, June 14-15, (2011).
%
\bibitem{Pfeifer07}
 A. A. Pfeifer, ``Controle Ótimo de Sistemas Algébrico-Diferenciais com Flutuação do Índice Diferencial'',
 Dissertação de Mestrado, Faculdade de Engenharia Química, Universidade Federal de Uberlândia, 2007.
%
\bibitem{Rashedi}
E. Rashedi, ``Gravitational Search Algorithm'',
M.Sc. Thesis, Shahid Bahonar University of Kerman, Kerman, Iran, 2007.
%
\bibitem{Schutz}
B. Schutz, {\em Gravity from the Ground Up}, Cambridge University Press, (2003).
%
\bibitem{SrPB03}
B. Srinivasan, S. Palanki and D. Bonvin,
Dynamic Optimization of Batch Processes: I - Characterization of the Nominal Solution, {\em Computers and Chemical Engineering}, {\bf 27}, 1-26, (2003).
%
\bibitem{Sakly09}
M. Sakly, A. Sakly, N. Majdoub and M. Benrejeb, Optimization of Switching Instants for Optimal Control of Linear Switched Systems based on Genetic
Algorithms, {\em 2nd IFAC International Conference on Intelligent control Systems and signal Processing}, 1-5, ICONS 2009, Istanbul, Turkey, (2009).
%
\bibitem{Sakly11}
M. Sakly, A. Sakly, N. Majdoub and M. Benrejeb, Optimal Switching Instants of Linear Switched Systems Based on Meta-heuristic Approaches,
{\em International Journal of Intelligent Control and Systems}, {\bf 16}(1), 8-18, (2011).
%
\bibitem{Xu01}
X. Xu, ``Analysis and Design of Switched Systems'', D.Sc. Thesis, University of Notre Dame, 2001.
%
\bibitem{Xu03}
X. Xu and P. J. Antsaklis, Results and Perspectives on Computational Methods for Optimal Control of Switched Systems. {\em Hybrid Systems: Computation and
Control}, Prague, Czech Republic, {\bf 2623}, 540-555, (2003).

\bibitem{XuAn04}
X. Xu and P. J. Antsaklis, Optimal Control of Switched Systems based on Parameterization of	the Switching Instants, {\em IEEE Transactions on Automatic Control}, {\bf 49}, 1-16, (2004).
\end{thebibliography}

\end{document}
\newpage
$ \  \  $  \thispagestyle{myheadings}  \markboth{      }{   }