ads:
ABORDAGEM DE UM PROBLEMA DE CUSTOS
PENALIZADOS NA DISTRIBUIÇÃO DE GÁS NATURAL
FORMULADO POR UM MODELO DE PROGRAMAÇÃO
MATEMÁTICA EM DOIS NÍVEIS
SARA MEIRA MOUTTA
UNIVERSIDADE ESTADUAL DO NORTE FLUMINENSE
DARCY RIBEIRO - UENF
CAMPOS DOS GOYTACAZES/RJ
FEVEREIRO – 2008
ads:
Livros Grátis
http://www.livrosgratis.com.br
Milhares de livros grátis para download.
II
Abordagem de um problema de Custos Penalizados na
distribuição de Gás Natural formulado por um modelo de
Programação Matemática em Dois Níveis
Sara Meira Moutta
"Dissertação apresentada ao Centro de
Ciência e Tecnologia, da Universidade
Estadual do Norte Fluminense, como parte
das exigências para obtenção do título de
Mestre em Engenharia de Produção, na
área de concentração de Pesquisa
Operacional.”
Orientadora: Profª. Gudelia G. Morales de Arica
Universidade Estadual Do Norte Fluminense Darcy Ribeiro – Uenf
Campos dos Goytacazes, RJ
Fevereiro, 2008
ads:
III
Abordagem de um problema de Custos Penalizados na
distribuição de Gás Natural formulado por um modelo de
Programação Matemática em Dois Níveis
Sara Meira Moutta
"Dissertação apresentada ao Centro de
Ciência e Tecnologia, da Universidade
Estadual do Norte Fluminense, como parte
das exigências para obtenção do título de
Mestre em Engenharia de Produção, na
área de concentração de Pesquisa
Operacional.”
Aprovada em 27 de fevereiro de 2008.
Comissão Examinadora:
______________________________________________
Dr. Edson Kenji Iamashita - PETROBRÁS
______________________________________________
Prof. José Ramón Arica Chávez – LEPROD/CCT/UENF
______________________________________________
Prof. Luis Humberto Guillermo Felipe – LCMAT/CCT/UENF
______________________________________________
Profª. Gudelia G. Morales de Arica – LEPROD/CCT/UENF
Orientadora
IV
Aos meus pais, Daniel e Eliene,
dedico esta obra.
V
AGRADECIMENTOS
A Deus por tudo o que alcancei até aqui, pois tenho certeza de que Ele sempre
esteve comigo, guiando-me a cada passo desta longa caminhada.
À minha orientadora, a professora Gudelia, pelo incentivo, paciência e todas as
horas dedicadas a me orientar na elaboração deste trabalho, principalmente os finais
de semana.
Ao professor Luis Guillermo, um grande amigo desde a graduação, pelo apoio e
ensinamentos transmitidos.
Aos autores do artigo que originou este trabalho, Dr. Kalashnikov e Dr. Ríos-
Mercado, pela atenção e boa vontade em compartilhar os seus conhecimentos.
Aos professores e funcionários da UENF, em especial o professor Arica e os
funcionários Kátia e Rogério, pela amizade e apoio.
Aos colegas do mestrado que me acompanharam nesta batalha, tais como o
Gilberto, o Luis Henrique, a Pollyanna e o Paulão. Sabemos que o caminho é árduo,
mas quando temos amigos ao nosso lado, tudo fica bem melhor.
Aos amigos da Fenorte que sempre me incentivaram neste mestrado, em especial a
Valéria Ferreira, que no momento em que precisei de liberação para estudar, ainda
como aluna especial, ela abraçou a minha causa e confiou em mim.
Aos amigos do Cederj, do C. E. Benta Pereira, e todos os demais amigos que
carinhosamente me apoiaram e acreditaram que seria possível.
Aos meus familiares, em especial os meus pais, que com muito amor e sabedoria
me ajudam a crescer e lutar pelos meus sonhos.
Ao meu noivo, Guilherme, por estar ao meu lado em todos os momentos.
VI
SUMÁRIO
Resumo ...........................................................................................................XVII
Abstract ..........................................................................................................XVIII
Capítulo …………………………………………………………………………………1
INTRODUÇÃO .................................................................................................... 1
Capítulo 2 ............................................................................................................ 6
UMA FORMULAÇÃO DO PROBLEMA DE DESBALANÇO DOS VOLUMES DE
GÁS NATURAL EM GASODUTOS ..................................................................... 6
2.1 – O Modelo Matemático do Problema ........................................................... 6
2.1.1 - Notação .................................................................................................... 8
2.1.2 - Problema de Nível Superior ..................................................................... 9
2.1.3 - Problema de Nível Inferior...................................................................... 10
2.1.4 - Comentários das restrições do 1° nível:................................................. 11
2.1.5 - Comentários das restrições do 2° nível:................................................. 11
Capítulo 3 .......................................................................................................... 14
FUNDAMENTAÇÃO TEÓRICA ......................................................................... 14
3.1 – Desigualdades Variacionais ..................................................................... 14
3.2 – Multiplicadores de Lagrange..................................................................... 16
3.2.1 – A função Lagrangeana na definição de dualidade de um problema de
otimização ......................................................................................................... 19
3.3 – Funções Gap ............................................................................................ 20
3.3.1 – Propriedades da Função Gap................................................................ 21
3.3.2 – Formulação do Ponto Sela .................................................................... 24
3.3.3 – Função Gap Primal e suas propriedades .............................................. 25
3.3.4- Função Gap Dual e suas propriedades................................................... 27
3.3.5 – Exemplo de uma função Gap para um programa Convexo................... 28
3.4 – Método de Penalização ............................................................................ 30
3.4.1 – O conceito de Função Penalizada ......................................................... 30
3.4.2 - Penalidade Exata e Inexata para modelos de Otimização com Restrições
.......................................................................................................................... 32
3.4.3 - Penalidade Exata e Inexata para problemas de otimização com uma
Desigualdade Variacional.................................................................................. 34
Capítulo 4 .......................................................................................................... 38
PRÉ-PROCESSAMENTO ................................................................................. 38
4.1 – Reformulação do Modelo do Problema de Dois Níveis ............................ 38
4.2 – Reformulação das restrições do problema de dois níveis ........................ 42
Capítulo 5 .......................................................................................................... 44
EXPERIÊNCIAS NUMÉRICAS.......................................................................... 44
5.1 – Primeiro Problema teste ........................................................................... 44
5.2 – Segundo Problema teste .......................................................................... 47
5.3 – Terceiro Problema teste ........................................................................... 49
5.4 – Quarto Problema teste.............................................................................. 51
5.4 – O Problema das redes de transporte de gás natural ................................ 54
Capítulo 6 .......................................................................................................... 60
CONCLUSÕES ................................................................................................. 60
7 – Referências Bibliográficas ........................................................................... 62
Anexo A............................................................................................................. 65
Apêndice A ........................................................................................................ 69
VII
RESUMO
Este trabalho estuda o modelo de um dos muitos problemas que surgem com a
transmissão e negociação de gás natural, apresentado por Kalashnikov e Ríos-
Mercado (2006). O modelo reflete as dificuldades dos emissores (responsáveis pela
extração do gás) para entregarem a quantidade de gás estipulada por contratos
feitos com os transportadores (responsáveis pelo transporte do gás através de
dutos). As dificuldades aparecem devido à natureza da própria indústria de gás, e
são chamadas de desequilíbrios.
Para evitar desequilíbrios muito grandes, emissores e transportadores adotam a
cobrança de multas ou taxas proposta pelo transportador, aceita por ambas as
partes. O emissor, por sua vez, busca tomar decisões de forma que minimize suas
multas, e, conseqüentemente, seus custos. O problema de minimizar estes custos é
tratado neste trabalho e representado por um modelo de programação matemática
em dois níveis.
Este trabalho busca entender, fundamentar e desenvolver uma alternativa
computacional de resolução do modelo proposto, reformulando o problema de
minimização do segundo nível como uma Desigualdade Variacional, cuja solução
encontrada através de um problema de otimização equivalente aproveita
propriedades da função de mérito associada, chamada função gap, facilitando a
solução do problema de dois níveis, via métodos de penalização.
Para verificar os métodos aplicados foram usados problemas testes e algoritmos de
resolução para estes problemas, implementados usando o ambiente MATLAB. Por
fim, apresentam-se os resultados de experiências computacionais, algumas
conclusões e trabalhos futuros a serem desenvolvidos sobre o problema abordado.
Palavras-chaves: programação em dois níveis, pré-processamento, desigualdade
variacional, função gap, método de penalização, distribuição de gás natural.
VIII
ABSTRACT
This work studies the model of one of the many problems that arise with the
transmission and marketing of natural gas: a natural gas cash-out problem, depicted
by Kalashnikov and Ríos-Mercado (2006). The model reflects the difficulties of the
shipper (responsible for extraction of gas) to deliver the exact amount of gas
contracted made to the pipeline company (responsible for transportation of the gas
through the pipes). The difficulties show up due to the nature of the gas natural
industry itself, and are called imbalance.
In order to avoid too big imbalance, shippers and pipeline companies adopt
the cash-out penalty policy, agreed by both parts. The shipper, on its turn, tries to
make decisions in order to minimize its penalties, and therefore its costs. The
problem related to minimizing the cash-out penalty is seen on this work and is
represented by a bilevel mathematical programming model.
This work tries to comprehend, ground and developing a computational
alternative to solve the proposed model by reformulating the problem of second level
minimization with a variational inequality, whose solution is found through an
equivalent optimization problem using properties of a merit function called gap
function which help to solve the bilevel programming problem with penalty method.
Numerical testing problems along with their resolution algorithms were used to
validate the applied methods. The algorithms were implemented using MATLAB.
Finally, the results from computational experiments were depicted, as long with some
conclusions and further works to be developed about the approached problem.
Keywords: bilevel programming, preprocessing, variational inequality, gap function,
penalty method, natural gas distribution.
1
Capítulo 1
INTRODUÇÃO
Este trabalho estuda um modelo de programação de dois níveis, que aborda um dos
muitos problemas que surgem com a transmissão e negociação de gás natural. O
problema surge no sistema de transporte de gás natural, onde se podem destacar
dois importantes subsistemas: os emissores que estão incumbidos da extração do
gás; e as empresas responsáveis pelo gasoduto, os transportadores, que estão
incumbidos do transporte do gás. Para realizar o transporte do gás, os
transportadores usam uma rede de tubulações que são chamadas de gasodutos. O
gás que flui através dos gasodutos perde energia devido à resistência das paredes
dos mesmos, e a pressão do gás vai caindo ao longo das tubulações. Sendo assim,
para que o gás chegue ao seu destino final, é preciso que se coloquem
compressores em alguns pontos dos gasodutos, para que a pressão necessária seja
mantida e, assim, o gás continue fluindo, como acontece com qualquer fluido em
movimento. Os pontos onde os compressores são colocados constituem o que é
chamado de estações de compressão. Estas são entidades complexas, subsistemas
menores, que envolvem um número de unidades de compressão organizadas em
paralelo ou em série, em uma única posição. Entretanto, para que os compressores
funcionem, há um gasto de combustível, pois estes consomem parte do gás natural
para se manterem em funcionamento. O custo operacional da rede de transporte de
gás é determinado, geralmente, pelo custo operacional das estações de
compressão, que varia de acordo com a pressão de sucção, a pressão de descarga
do gás e o fluxo que passa pelo duto. Sendo assim, o gás a ser transportado gera
custos e, quanto maior o volume de gás a ser transportado, maior será o custo.
O problema de minimizar estes custos é motivo de estudo de muitos pesquisadores.
Busca-se a sua solução através de diferentes caminhos. Uma das estratégias
utilizadas é o uso de técnicas conhecidas para manipular o modelo do problema,
tornando-o mais apropriado para os métodos de solução existentes, mantendo
rigorosamente a equivalência com o modelo original, o que chamamos de pré-
processamento.
2
Wu et al. (2000) apresentaram um método de busca de solução do problema,
através de um esquema de limitação inferior, baseado na técnica de relaxação
trabalhados nas seguintes situações: relaxação da função objetivo do custo do
combustível usado pelo compressor e a relaxação do conjunto viável (o domínio);
ambas originalmente são não-lineares e não-convexas.
Segundo Rios-Mercado et al. (2000), o problema foi modelado de modo que a
função objetivo seja o custo total de combustível utilizado nas estações
compressoras da rede. As funções envolvidas no modelo são tipicamente não-
convexas, não-lineares e não diferenciáveis. Eles resolveram utilizar a metodologia
resultante da combinação da teoria dos grafos com análise funcional não-linear e
conseguem reduzir significantemente a dimensão das variáveis e/ou o número de
restrições do problema em questão.
Ríos-Mercado et al. (2006) fizeram uso de uma metodologia de cálculo de uma
solução heurística, baseada em um procedimento iterativo envolvendo dois estágios.
No primeiro, as variáveis do fluxo do gás são fixas e as variáveis ótimas de pressão
são encontradas através da programação dinâmica. No segundo, as variáveis de
pressão são fixas e é feita uma tentativa para encontrar um conjunto das variáveis
do fluxo, que aprimoram a função objetivo através da exploração da estrutura básica
da rede.
Kalashnikov e Ríos-Mercado (2006) formularam um modelo que será tratado neste
trabalho. O modelo reflete as dificuldades de um emissor para entregar a exata
quantidade de gás contratada. Inicia-se o problema quando um emissor faz um
contrato com uma companhia de gasodutos (transportador) para entregar uma
determinada quantidade de gás em diversos pontos de entrega. Devido à natureza
da própria indústria de gás, a quantidade de gás entregue na prática pode ser maior
ou menor do que o volume contratado, o que acarreta um desequilíbrio. Para evitar
desequilíbrios muito grandes, emissores e transportadores adotam a cobrança de
multas ou taxas proposta pelo transportador, aceita por ambas as partes, para
quando for o caso de haver um desequilíbrio. O emissor, por sua vez, busca tomar
decisões de forma que minimize suas multas. No ambiente de comercialização do
3
gás natural, os responsáveis pelas negociações antes mencionadas, descrevem
este mecanismo como cash-out penalty.
O procedimento de tomada de decisão dos emissores e transportadores de gás
natural é interdependente, onde a tomada de decisão de cada um destes influencia
o comportamento e a tomada de decisão do outro.
No processo de tomada de decisão descrito anteriormente, existe uma hierarquia
responsável pela tomada de decisão, que é feita em níveis diferentes e
complementares. Um modo para representar tais hierarquias é formulando níveis,
isto é, para cada hierarquia deve-se definir um nível e uma função de decisão e suas
respectivas variáveis e restrições. Quando as funções envolvidas em todos os níveis
são lineares, tem-se um modelo de programação matemática linear em níveis.
Apresenta-se a seguir o modelo genérico de programação matemática de dois
níveis, que é o modelo que analisa problemas de tomada de decisão considerando
duas estruturas hierárquicas.
A formulação geral de um modelo Dois Níveis é:
Min. F(x, y)
s.a.
g
i
(x, y) ≤ b
i
, i = 1, 2, ..., k (1.1)
e y solução de
Min. f(x, y)
s.a.
h
j
( x,y) ≤ d
j
, j = 1, …, p (1.2)
(x, y) ∈ D
Onde as funções F(x, y), f(x, y), g
i
(x, y) e h
j
(x, y) são funções definidas nas
dimensões apropriadas e D é o conjunto de definição do problema.
O conjunto das soluções viáveis associado com um problema de programação de
dois níveis, chamado conjunto induzido na literatura, é determinado
4
seqüencialmente por dois problemas de otimização, que devem ser resolvidos em
estágios predeterminados:
- Temos um líder (associado com o nível superior) que seleciona sua decisão
primeiramente, associado a variável x;
- E o seguidor (associado com o nível inferior) que responde a esta decisão.
Em outras palavras, as incógnitas deste modelo (x, y), são divididas em duas
parcelas onde x é considerada a variável do líder e y do seguidor.
Usando x como um parâmetro, o seguidor resolve um problema paramétrico de
otimização que computa uma solução ótima y*= y*(x). Esta solução ótima y*(x) é a
reação (ou resposta) do seguidor para a seleção x do líder. Sabendo as reações do
seguidor em suas decisões, o líder tem agora que minimizar uma determinada
função F(x,y*(x)) a respeito de suas variáveis de decisão.
O problema de transporte de gás natural citado acima é representado por um
modelo particular de programação matemática em dois níveis. Ele foi modelado por
Kalashnikov e Ríos-Mercado (2006), identificando funções lineares com variáveis
inteiras e contínuas, constituindo-se em um modelo de dois níveis misto-inteiro, onde
o emissor será representado pelo líder e o transportador será representado pelo
seguidor.
O estudo e compreensão desta proposta de modelo contribuem para que se tenha
um conhecimento mais amplo de uma metodologia e fundamentos teóricos e
computacionais, pelos quais o problema das redes de transporte de gás natural
pode ser abordado e tratado com o objetivo da busca de soluções ótimas.
A proposta deste trabalho é entender, fundamentar e desenvolver uma alternativa
computacional de resolução do modelo proposto, reformulando o problema de
minimização do segundo nível como uma desigualdade variacional, cuja solução
será encontrada através de um problema de otimização equivalente, que usa uma
função de mérito chamada função gap. Toda esta manipulação do modelo faz com
que o problema de programação de dois níveis se transforme em um problema de
programação matemática padrão, sem que suas restrições sejam alteradas. Para
5
encontrar a solução deste modelo equivalente, já existem algoritmos testados e
ambientes computacionais confiáveis, o que facilitará a resolução do modelo do
problema.
A estrutura do trabalho é a seguinte: no Capítulo 2 apresenta-se a formulação do
problema e o modelo apresentado por Kalashnikov e Ríos-Mercado (2006). No
Capítulo 3 apresenta-se uma fundamentação teórica onde são apresentadas: a
desigualdade variacional, as potencialidades dos multiplicadores de Lagrange, fonte
de inspiração para a definição da função gap, e o método de penalidades para o
cálculo de uma solução. No Capítulo 4 apresenta-se o pré-processamento aplicado
no modelo formulado no capítulo 2, com o objetivo de facilitar o cálculo de solução, e
algumas reformulações deste modelo feitas visando adequá-lo ao ambiente
MATLAB. No Capítulo 5 apresentam-se os resultados da experiência
computacional, da implementação de um algoritmo em ambiente MATLAB. E,
finalmente, no Capítulo 6 apresentam-se algumas conclusões e futuros trabalhos a
serem desenvolvidos sobre o problema abordado.
6
Capítulo 2
UMA FORMULAÇÃO DO PROBLEMA DE DESBALANÇO DOS VOLUMES
DE GÁS NATURAL EM GASODUTOS
Neste capítulo apresenta-se uma formulação do problema de minimizar os custos
por desequilíbrios contratuais de entrega de gás natural, valendo-se de multas sobre
as diferenças contratuais dos volumes efetivos de gás (cash-out prices). Utilizou-se a
proposta de Kalashnikov e Ríos-Mercado (2006) que modela o problema usando
programação matemática em dois níveis.
2.1 – O Modelo Matemático do Problema
No sistema de transporte de gás natural existem os emissores (responsáveis pela
extração do gás) e as companhias de gasodutos ou transportadores (responsáveis
pelo transporte do gás). Estes assinam contratos que estipulam a quantidade de gás
a ser transportada pelos gasodutos. Para que se tenha o controle desta quantidade,
são instalados medidores de entrega nos pontos onde os emissores colocam o gás
a ser entregue e também medidores de recebimento nos pontos onde os
transportadores recebem este gás.
Suponha que um emissor assinou um contrato para entregar uma certa quantidade
de gás natural a ser mensurado por um medidor de recebimento em um determinado
tempo. O emissor estipula certa quantidade diária de gás a ser injetada pelo
operador do medidor de entrega e a ser retirada pelo operador do medidor de
recebimento. O gasoduto transporta o gás do medidor de entrega ao medidor de
recebimento. Devido à natureza da indústria do gás natural, o que é realmente
transportado é inevitavelmente diferente do que é contratado. Tal diferença constitui
um desequilíbrio ou desbalanceamento.
Existem desequilíbrios operacionais e de transporte. O primeiro tipo de desequilíbrio
refere-se às diferenças entre fluxos nominais e os fluxos reais. O desequilíbrio de
7
transporte envolve diferenças entre a quantidade do gás inicialmente entregue,
colocado no medidor de entrega, e a quantidade do gás ao finalizar o transporte,
aquele que sai no medidor de recebimento. Neste caso, a diferença que se pode
observar, vem do gasto de gás como combustível para transitar de um medidor a
outro.
Quando acontecem pequenos desequilíbrios, não há com o que se preocupar,
porém quando eles são consideráveis, são emitidas taxas ou multas. Estas são
criadas em comum acordo entre emissores e transportadores com a finalidade de
evitar desbalanceamentos muito grandes. Estas taxas consistem no pagamento de
multas por parte do emissor, para o caso de haver um desequilíbrio que ultrapasse
os limites estipulados neste acordo. Assume-se, no modelo, que os emissores foram
informados e são cientes das taxas de desequilíbrio estabelecidas no contrato e se
analisará o pagamento de taxas associado com as operações de desequilíbrio.
Por parte do emissor, um desequilíbrio operacional pode ser positivo ou negativo.
Um desequilíbrio positivo acontece quando o transportador retira mais gás do que o
combinado. Um desequilíbrio negativo acontece quando o transportador retira
menos gás do que o combinado. Desequilíbrios positivos [negativos] no fim de mês
implicam o pagamento da diferença efetuado pelo transportador [emissor]. As taxas
por desequilíbrio (cash-out prices) são ajustadas de uma maneira que sempre que o
emissor vende [compra] o gás para o [do] transportador, ele o faça a um preço baixo
[elevado]. A relação entre as taxas por desequilíbrio (cash-out prices) e a posição do
desequilíbrio é não linear com respeito ao valor médio entre os preços máximo e
mínimo do mesmo dia do mês anterior.
Assim como é verdade que os desequilíbrios se devem a fatores externos, também é
verdade que o emissor tem um importante papel nas decisões quanto aos
desequilíbrios. Pelo fato de que as componentes externas, caso sejam modeladas,
gerariam variáveis aleatórias, o que resultaria além do foco deste trabalho, foi
assumido que as decisões de desequilíbrio se devem totalmente ao emissor.
Portanto, modelou-se como o emissor deve atentar para os desequilíbrios e seus
custos. O preço de gás é um dos fatores principais que afetam as decisões dos
emissores. Na ausência de um acordo ou provisões de venda, historicamente o
emissor removeria o gás do gasoduto no inverno sob altos custos (que causa
8
desequilíbrios negativos), e pagariam o transporte de volta com custos baixos no
verão. Isto corresponde a um comportamento especulativo pelo emissor, por meio
do qual os desequilíbrios são criados e controlados a fim de tirar vantagem com os
movimentos do preço de gás. As taxas de venda foram projetadas a fim de evitar tais
arbitrariedades, fixando o preço.
2.1.1 - Notação
As notações seguintes são usadas para descrever o modelo:
Índices e conjuntos de índices
i,j,k índices de zona; i, j, k ∈ J = {1, 2,...,P};
t períodos do tempo; t ∈T = {1, 2,..., N}, sendo T o horizonte em estudo.
Parâmetros
x
ti
L
, x
ti
U
limites, inferior e superior respectivamente, dos desequilíbrios diários
no fim do dia t na zona i; t ∈T, i ∈J;
x
t
L
, x
t
U
limites totais inferior e superior dos desequilíbrios diários, no fim do dia
t; t ∈T;
s
ti
L
, s
ti
U
limites, inferior e superior, de variação de desequilíbrio durante o dia t
na zona i; t ∈T, i ∈J;
e
ij
porcentagem do combustível retido para movimentar um dt de gás da zona
i para j; i, j ∈ J. O termo dt é uma abreviação da palavra, de origem inglesa,
dekatherm. Significa uma energia calórica equivalente a 1 milhão de BTUs;
f
ij
custo do transporte (empurrar) de um dt de gás da zona i para j; i, j ∈ J, i < j;
b
ij
bonificação do transporte (puxar) de um dt de gás da zona j para i; i, j ∈ J, i
< j;
x
0j
desequilíbrio inicial (o modelo começa no dia 1) na zona j; j ∈J;
r
j
, δ
j
parâmetros de penalização pelo desequilíbrio na zona j ∈J.
Variáveis de decisão (primeiro nível)
9
x
ti
desequilíbrio no fim do dia t na zona i; t ∈T, i ∈J;
s
ti
variação do desequilíbrio durante o dia t na zona i; t ∈T, i ∈J.
Variáveis de decisão (segundo nível)
y
i
desequilíbrio final na zona i; i ∈J;
u
ij
volume transportado (empurrado, despachado) da zona i para j; i, j ∈ J, i <
j;
v
ij
volume transportado (puxado, devolvido) da zona j para i; i, j ∈ J, i < j;
z receita total, sob contrato, do emissor.
Variáveis auxiliares
q variável binária igual a 1 (0) se os desequilíbrios finais são não-negativos
(não-positivos). No caso especial, quando todos os desequilíbrios finais são
zeros, será aceito q = 1.
2.1.2 - Problema de Nível Superior
Este problema tem como objetivo minimizar a receita total contratual do emissor.
Esta receita é representada pela variável z, que será descrita em (2.2i).
Minimizar h
1
(x,s,y,u,v,z) = z (2.1a)
(x,s,y,u,v)
Sujeito a:
(2.1b)
(2.1c)
(2.1d)
(2.1e)
JiTtxxx
U
titi
L
ti
∈∈≤≤ ,
JiTtsss
U
titi
L
ti
∈∈≤≤ ,
Ttxxx
U
t
Ji
ti
L
t
∈≤≤
∑
∈
JiTtsxx
tiitti
∈∈+=
−
,
,1
10
Onde o vetor (y,u,v) é a solução ótima do problema de nível inferior, fixadas as
variáveis x e s, e está constituído por
y
, variável que representa os desequilíbrios
finais de gas em cada zona;
u,
variável de volume de despacho de gás entre zonas
e
v
, variável de volume de devolução de gás entre zonas.
2.1.3 - Problema de Nível Inferior
Sabe-se que as taxas cobradas quando aparecem desequilíbrios são um acordo
estabelecido por ambas as partes, emissor e transportador, e quanto maiores forem
estas taxas, maior será o custo do emissor. O objetivo deste nível não pode ser
maximizar a receita do emissor, pois assim seria um incentivo apenas para o
transportador. Neste caso, eles combinam em minimizar os desvios de zero da
função custo por desequilíbrios, isto é, minimizar a quantidade de transações de
capital positiva ou negativamente. Deste modo, o emissor ficará protegido contra
valores negativos da receita total do emissor e o transportador ficará protegido
contra valores positivos desta receita que é formulada na restrição (2.2i) abaixo.
Minimizar h
2
(x,s,y,u,v,z) = |z|, (2.2a)
(y,u,v)
sujeito a:
(2.2b)
(2.2c)
(2.2d)
(2.2e)
(2.2f)
<>
≤
ontrário. caso , 0
;0 xe 0 xse ,
jN,iN,
c
x
u
Ni
ij
<>
≤
ontrário. caso , 0
;0 xe 0 xse ,
iN,jN,,
c
x
v
jN
ij
∑ ∑∑∑
> <><
∈−−+−+=
jkk jii
ijjk
jkk
jk
jii
ijijjNj
Jjvuvuexy
: :::
,
;,)1(
∑
∑
∑
∈ < ∈
=+
Ji jiji Ji
iNijiji
xuey
:),(
,
Jixvu
iN
ijj ikk
kiij
∈≤+
∑ ∑
< <
},,0max{
,
: :
11
(2.2g)
sendo M um inteiro positivo. (2.2h)
(2.2i)
Onde (y
i
)
+
= max { 0,y
i
}.
y
i
, z variáveis sem restrição de sinal, i ∈ J; (2.2j)
(2.2k)
(2.2l)
2.1.4 - Comentários das restrições do 1°
°°
° nível:
As restrições formuladas em (2.1b) estipulam limites dos desequilíbrios diários
ao final do dia i , na zona j; visto que desequilíbrios muito grandes não
satisfazem as expectativas dos emissores nem dos transportadores;
As restrições em (2.1c) estabelecem os limites de variação do desequilíbrio
durante o dia t, pois deverá existir um controle das alterações do desequilíbrio
na zona j;
As restrições em (2.1d) asseguram que a soma dos desequilíbrios diários em
todas as zonas em um determinado dia deve estar entre os limites estipulados
para este dia;
As restrições em (2.1e) garantem que em qualquer zona i, o desequilíbrio ao
final de um dia t deve ser igual ao desequilíbrio ao final do dia anterior,
somado com a variação de equilíbrio durante o dia.
2.1.5 - Comentários das restrições do 2°
°°
° nível:
As restrições formuladas em (2.2b) representam a relação entre o
desequilíbrio ao final do dia N na zona j, os volumes de gás puxado ou
{
}
{
}
.,,0max,0min
,,
Jixyx
iNiiN
∈≤≤
,,)1( JiMqyqM
i
∈≤≤−−
[
]
∑
∑
∑
∈ < <
+
−+−−=
Ji jiji jiji
ijijijijijiiii
uefvbyryz
:),( :),(
.
2
)1()(
δ
;0, ≥
ijij
vu
{
}
.1,0
∈
q
12
empurrado além dos volumes consumidos de gás que interagem com a zona j
e o desequilíbrio final na zona j;
As restrições em (2.2c) asseguram que nenhuma perda de gás ocorrerá;
As restrições em (2.2d) decidem sobre o que será movimentado em cada
zona, impedindo movimentos cíclicos de gás. Isto simplesmente afirma que,
a partir de qualquer zona dada i, não se pode mover mais do que o
desequilíbrio positivo inicial;
As restrições em (2.2e) e (2.2f) asseguram o sentido certo de fluxo do gás,
empurrado e/ou puxado entre as zonas i e j, relacionado aos desequilíbrios
das respectivas zonas ao final do último dia, N;
As restrições em (2.2g) asseguram que os desequilíbrios finais na zona i, y
i
, e
os desequilíbrios no último dia no final do dia da zona i,
x
N,i
, terão o sinal
certo, isto é, ambos devem ter o mesmo sinal;
As restrições em (2.2h) mostram que os desequilíbrios finais para todas as
zonas devem ter o mesmo sinal, ou seja, os desequilíbrios não devem mudar
de sinal;
A restrição em (2.2i) representa o custo do ponto de vista do emissor. Dividiu-
se z em três parcelas. Na primeira se incluem os parâmetros de multa ou
taxas contratuais para os desequilíbrios totais por zona; na segunda se
representa o crédito do transporte de gás “devolvido” ou puxado da zona j
para a zona i; e na terceira parcela representa o gasto com o transporte de
gás enviado ou empurrado da zona i para a zona j. Observe que se y
i
> 0, isto
é, quando ocorre um desequilíbrio positivo, a restrição (2.2i) torna-se:
(2.2i a)
É possível ver que a primeira parcela de (2.2i a) não tem sinal definido,
podendo ser positiva ou negativa, de acordo com o valor do desequilíbrio total
na zona i, y
i
.
Por outro lado, se y
i
< 0, isto é, quando ocorre um desequilíbrio negativo, a
restrição (2.2i) torna-se:
(2.2i b)
Neste caso pode-se ver que a primeira parcela de (2.2i b) é positiva.
[
]
∑
∑
∑
∈ < <
−+−−=
Ji jiji jiji
ijijijijijiiii
uefvbyryz
:),( :),(
.
2
)1(
δ
[
]
∑
∑
∑
∈ < <
−+−−=
Ji jiji jiji
ijijijijijii
uefvbyrz
:),( :),(
.
)1(
13
As restrições em (2.2j) indicam que as variáveis y
i
e z podem assumir valores
positivos ou negativos, isto é, elas não têm restrição de sinal;
As restrições em (2.2k) indicam que as variáveis u
ij
e v
ij
não podem assumir
valores negativos;
A restrição em (2.2l) estabelece que a variável q só pode assumir os valores:
q=0 ou q=1.
Maiores comentários sobre as restrições do 2° nível podem ser vistos em
Kalashnikov e Ríos-Mercado (2006).
14
Capítulo 3
FUNDAMENTAÇAO TEÓRICA
Neste capítulo serão apresentados conteúdos matemáticos necessários para a
compreensão do modelo formulado e o entendimento da metodologia computacional
utilizada no cálculo de soluções do modelo. Dentre os conteúdos estudados, podem-
se destacar as Desigualdades Variacionais, os Multiplicadores de Lagrange, as
funções Gap e o método de Penalidades. Estes conteúdos, apresentados de
maneira mais detalhada nas seções seguintes, formam toda a base teórica na qual
se apóia a resolução de um modelo equivalente ao modelo apresentado no capitulo
2, que é numericamente menos custoso, cujos detalhes serão apresentados no
capítulo 4.
3.1 – Desigualdades Variacionais
Historicamente, Hartman e Stampacchia (1966) (apud Trujillo (2003)) introduziram a
teoria das desigualdades variacionais como uma ferramenta para o estudo de
problemas mecânicos formulados usando equações diferenciais parciais. Entretanto,
um dos eventos que marcou o início da teoria das desigualdades variacionais
aconteceu na década de 80, quando Dafermos (1980) reconheceu que uma
condição de equilíbrio em redes de transporte, formulada por Smith (1979),
apresentava a estrutura de uma desigualdade variacional. Daí em diante a teoria das
desigualdades variacionais foi amplamente utilizada para o estudo das condições de
equilíbrio em procedentes de diferentes problemas, tais como: da economia, da
ciência administrativa/pesquisa operacional, da física, e da engenharia do tráfego.
Definição 3.1: Seja um conjunto C ⊆ R
n
,
não vazio, convexo e fechado, e F: C → R
n
uma função contínua dada. O problema de desigualdade variacional de F sobre C,
PDV(F,C), é definido como:
Encontrar um vetor x
∈
C que satisfaz
〈
F(x), (y - x)
〉
≥ 0,
∀
y
∈
C (3.1)
15
Observar que sendo C ⊆ R
n
e F: C → R
n
Se n = 1, a desigualdade variacional torna-se:
F(x)(y-x) ≥ 0.
Se n = 2, F = (F
1
, F
2
), a desigualdade variacional torna-se:
〈
(F
1
(x), F
2
(x)),((y
1
– x
1
), (y
2
– x
2
))
〉
≥ 0 que é equivalente a
F
1
(x)(y
1
– x
1
) + F
2
(x)(y
2
– x
2
) ≥ 0
Se n ≥3, a desigualdade variacional se escreve, então:
∑
=
≥−
n
i
iii
xyxF
1
0))((
(3.2)
Geometricamente, a relação (3.1) diz: deve-se ser encontrado um vetor x
∈
C
tal que o vetor imagem, F(x), faça um ângulo agudo com o vetor (y – x),
formado para cada y ∈ C.
Existem vários métodos para resolver um problema de desigualdade variacional. Um
deles obtém-se o reformulando para um problema de otimização equivalente (neste
caso, minimização), de apenas um nível, através de uma função associada. A
função objetivo do modelo de otimização associada ao PDV(F,C) é um indicador do
decrescimento (crescimento) da proximidade dos pontos candidatos ao ponto
solução do PDV(F, C). Assim, vai ser chamada neste trabalho de função de mérito.
A função de mérito associada ao PDV(F,C), formulada em (3.1), foi chamada de
função gap (Hearn (1981), Auchmuty (1989), Morales (1997)) e será estudada mais
detalhadamente na seção 3.3 deste capitulo.
Um caso particular do problema de desigualdade variacional é estabelecido, quando
a função F é o gradiente de alguma funcional f: R
n
→ R, convexa e diferenciável, isto
é, F(x) = ∇ f(x). Neste caso, o problema pode ser escrito como:
Encontrar x ∈ C, tal que, 〈∇ f(x),(y - x) 〉 ≥ 0, ∀ y ∈ C. (3.3)
Notar que a relação (3.3), sob as hipóteses acima indicadas para f, não é mais do
que a condição necessária e suficiente de otimalidade de primeira ordem.
16
3.2 – Multiplicadores de Lagrange
Os estudantes de Engenharia e Ciências: Física, Química ou Matemática são
apresentados aos multiplicadores de Lagrange como uma ferramenta para o cálculo
de soluções para modelos de problemas de otimização com restrições.
Originalmente foram introduzidos como um artifício algébrico na busca de uma
solução para modelos de otimização com restrições de igualdade, convertendo-se
em tema de pesquisa da otimização não-linear na década de 60, conseguindo-se
uma grande variedade de importantes resultados teóricos que ainda não terminaram
de ser explorados computacionalmente.
Inicia-se a apresentação do papel dos multiplicadores de Lagrange, considerando o
modelo de minimização com restrições de igualdade seguinte:
minimizar f(x)
sujeito a: h(x) = 0 (3.4)
onde as funções f : R
n
→R e h : R
n
→R são contínuas e diferenciáveis.
Geometricamente, isso significa que estamos olhando para um ponto x
0
, sobre o
gráfico da curva de restrição h(x
0
) = 0, na qual f(x) é tão pequeno quanto possível.
Para ajudar a localizar tal ponto, basta fazer algumas curvas de nível de f(x) = c.
Cada ponto de intersecção do lugar geométrico definido por h(x) = 0 com uma curva
de nível de f é um candidato para uma solução, uma vez que esses pontos situam-
se na curva de restrição. Entretanto, descobriu-se que o valor ótimo, neste caso o
mínimo de f(x), ocorre no ponto onde a curva de restrição e a curva de nível somente
se tocam, e neste ponto de tangência estas curvas têm uma reta normal em comum.
Isso sugere que o mínimo de f(x), se existir, ocorre em um ponto (x
0
) sobre a curva
de restrição na qual os vetores gradientes ∇f e ∇h são múltiplos escalares um do
outro, isto é,
∇f(x) = λ∇h(x) (3.5)
Para algum escalar λ. Este escalar é chamado de multiplicador de Lagrange.
Aumentando o número de restrições de igualdade e reescrevendo (3.4) tal que, as
17
funções h
i
são diferenciáveis numa vizinhança de x*, isto é, h
i
continuamente
diferenciável em x*, temos:
minimizar f(x)
sujeito a: h
i
(x) = 0, i = 1, . . . , m.
(3.6)
Adicionalmente é necessário garantir ou que o conjunto {∇h
i
(x*), i= 1, 2, 3, ...,m} seja
de vetores linearmente independentes ou que a função h seja afim linear (isto é,
h(x) = Ax - b com A∈R
n×m
e b∈R
m
) para garantir a extensão da relação (3.5). As
condições mencionadas, no parágrafo anterior, foram estudadas nos anos 60 e são
conhecidas como condições de qualificação.
O teorema de multiplicadores de Lagrange básico para o modelo (3.6) estabelece
que em um ponto mínimo local x* se verifica a relação:
∑
=
=∇+∇
m
i
ii
xhxf
1
0*)(*)(
λ
(3.7)
Da relação (3.7) se deduz que o vetor gradiente da função f, em um ponto de mínimo
local do modelo (3.6) é ortogonal ao subespaço de variações viáveis de primeira
ordem, representado pelo conjunto V (x*) = {dx | (∇h
i
(x*)dx) = 0; i = 1, . . . , m}. Esta
última observação pode ser estabelecida na seguinte proposição, Lema 2.2.1 de
Izmailov e Solodov (2005), que é um resultado básico da Álgebra Linear, reescrita a
seguir.
Proposição 3.1 -
Para qualquer matriz
A
∈R
n×m
se verifica que
(núcleo(
A
))* = imagem(
A
T
) = (núcleo(
A
))
⊥
.
A notação (K)* significa cone dual do cone K definido por K*= {y ∈ R
n
/ 〈y,k〉 ≤ 0,∀ k
∈K}; sendo núcleo(
A
) = {z
∈
R
n
|
A
z = 0}. No texto Izmailov e Solodov (2005) se usa
a notação ker(
A
) para descrevê-lo, im(
A
T
) para o conjunto das imagens resultantes
de aplicar a matriz
A
T
e (
S
)
⊥
indica o subespaços ortogonal ao subespaço
S
em R
n
.
O modelo (3.6) sob as condições de qualificação acima podem ser escrito como:
18
-∇f(x*) ∈ {núcleo[∇h(x*)] }* = imagem[∇h(x*)]
T
= {núcleo[∇h(x*)]}
⊥
,
onde [∇
h
(x*)] é a matriz jacobiana de
h
, o sistema de equações, restrição do
problema de minimização definido na relação (3.6). Pela proposição 3.1, pode-se ter
a visão geométrica de encontrarmos o negativo de ∇f(x*) no espaço ortogonal
definido pelos gradientes das restrições h
j
.
Este resultado é geralmente apresentado em termos do Lagrangeano do modelo
(3.6), igualado a zero. Sendo o Lagrangeano definido como uma função:
L: R
n
×
R
m
→R, cuja regra de correspondência é L(x, λ) =
∑
=
+
m
1i
ii
xhxf )()(
λ
. Então, a
condição (3.7) formula um sistema de n+m incógnitas e n+m equações, os valores
da variável x e dos multiplicadores λ a saber:
0xhxfxL
i
m
1i
ix
=∇+∇=∇
∑
=
)()(),(
λλ
e
∇
λ
L(x,
λ
) =h(x) = 0.
Também, ao considerar o modelo de otimização com restrições de desigualdade
minimizar f(x)
sujeito a: g
j
(x) ≤ 0, j = 1, . . . ,k. (3.8)
os multiplicadores de Lagrange podem ser visualizados como os preços de equilíbrio
do problema de otimização associado, o problema dual, sob as hipóteses das
funções f e g
j
, j=1,...k diferenciáveis e convexas sobre R
n
e o valor ótimo finito. Gale
(1967) (apud Ozdalgar (2003)) afirma ser esta interpretação a mais importante
ferramenta na análise econômica moderna, desde o ponto de vista teórico e
computacional simultaneamente.
Para o modelo (3.8), a condição clássica dos multiplicadores de Lagrange em um
mínimo global x*, sob condições apropriadas (veja, por exemplo, as condições de
qualificação na pagina anterior), diz que existem multiplicadores não negativos µ’
1
,
µ’
2
, ...,µ’
k
tal que
0*)(*)(),(
1
'
=∇+∇=∇
∑
=
xgxfxL
j
k
i
j
x
µλ
(3.9)
e os multiplicadores µ’
i
satisfazendo a condição de complementaridade de folgas
µ’
j
g
j
(x*) = 0. (3.10)
É conveniente lembrar que a hipótese de convexidade das funções envolvidas no
modelo com restrições de desigualdade torna a condição dos multiplicadores de
Lagrange em uma condição suficiente para o vetor mínimo global x*, da função
19
f(x)+
)(
1
,
xg
j
k
j
j
∑
=
µ
, isto é
,
com a condição de complementaridade de folgas para os
multiplicadores µ’
j
, o ponto x* é um ponto ótimo:
+=
∑
=
ℜ∈
k
j
jj
u
xgxfxf
n
1
)(')(inf*)(
µ
(3.11)
3.2.1 – A função Lagrangeana na definição de dualidade de um problema
de otimização
Consideramos agora um problema de minimização em formato geral (Izmailov e
Solodov, 2005) onde se apresentam restrições de igualdade e desigualdade.
minimizar f(x)
sujeito a: g
j
(x) ≤ 0, j = 1, . . . ,k.
h
i
(x) = 0, i = 1, . . . ,m. (3.12)
onde Ω ⊂ R
n
é um conjunto aberto e f : Ω → R, h : Ω → R
m
e g : Ω → R
k
são funções
convexas e inferiormente semi-contínuas (ver anexo A.3).
Seja D = {x ∈ Ω / h(x) = 0, g(x) ≤ 0} o conjunto viável do problema (3.12). Este
problema será chamado de problema primal.
A função Lagrangiana deste problema é dada por L : Ω x R
k
x R
m
→ R,
)()()(),,(
11
xgxhxfxL
j
k
j
ji
m
i
i
∑
∑
==
++=
µλµλ
(3.13)
Onde λ∈R
m
e µ∈R
+
k
A relação (3.13) pode ser reescrita usando o produto interno:
〉
〈
+
〉
〈
+
=
)()(xfxL xgµ,xh,λ)(),,(
µ
λ
. (3.14)
Se x ∈ Ω é viável, isto é, x ∈ D, temos que h(x) = 0 e g(x) ≤ 0. Logo se conclui que:
Sup L (x, λ, µ) = (3.15)
(λ, µ)∈ R
m
x R
k
+
Pode-se substituir (3.12) por:
sujeito a x ∈ D’. (3.16)
Ω∈∞+
∈
Dxse
Dxsexf
\,
),(
+
∈
∈
),,(supmin
),(
'
µλ
µλ
xL
km
xRR
Dx
20
onde D’
O problema dual é definido por:
)},,(inf{max
),(
µλ
µλ
xL
x Ω∈
∆∈
(3.17)
Onde ∆ = {(λ,µ) ∈ R
m
x R
k
+
/ inf L (x, λ, µ) > - ∞ }
Ao representar a função objetivo do problema dual, pela função valor
),,(inf),(
µλµλϕ
xL
x Ω∈
= , isto é, ϕ:∆→R;
o problema dual pode ser escrito como:
(3.18)
Esta formulação indica que no ponto (λ*, µ*)∈ R
m
×R
+
k
, mais as condições de
complementariedade, que relaciona o vetor
µ
com as restrições de desigualdade em
x*, fornecerão o equilíbrio entre os custos do modelo primal (3.12) e os lucros de seu
modelo dual.
3.3 – Funções Gap
Como um método alternativo para resolver problemas de desigualdades
variacionais Auslender (1976) (apud Trujillo (2003)), reformulou o problema de
desigualdade variacional PDV (F, C) apresentado em (3.1) como um problema
equivalente de otimização. Para isto, foi utilizada uma função de mérito, definida por:
〉−〈=
∈
yxxFxG
Cy
),(sup)( (3.19)
Onde o conjunto C ⊆ R
n
é não vazio, convexo e fechado, e F: C → R
n
é uma função
vetorial contínua dada. A função G(x) definida acima é chamada de função gap.
Se o conjunto C adicionalmente for limitado, pode-se reformular G(x) como:
(3.20)
No caso particular onde F(x) = ∇ f(x), para alguma f: R
n
→ R convexa e diferenciável,
a função gap associada ao problema de desigualdade variacional PDV (F, C) em
(3.3) é definida por:
)(),(min)( xyxFxG
Cy
−−=
∈
)(),(max)( yxxFxG
Cy
−=
∈
+∞<Ω∈=
+
∈
),,(sup
),(
µλ
µλ
xLx
km
xRR
),(max
),(
µ
λ
ϕ
µλ
∆∈
21
(3.21)
Em seguida será mostrada a propriedade mais útil, do ponto de vista computacional,
da função G(x) que é a de ser não negativa no conjunto C e ter o valor G(x) = 0 se, e
somente se, y∈C resolve o problema PDV (F, C) definido em (3.1). Assim, a
desigualdade variacional PDV (F, C) pode ser reformulada como um problema de
otimização equivalente onde será procurado o valor mínimo igual a zero:
Observação: na verdade, a minimização para a função gap em 3.21 é um problema
min
x
∈
C
max
y
∈
C
∇ f(x)
T
(x - y).
3.3.1 – Propriedades da Função Gap
Seja C ⊆ R
n
um conjunto não vazio, convexo e fechado, e F: C → R
n
uma função
contínua dada.
Considere a função Φ
:
C → R ∪ {∞} definida, Auchmuty (1989), por:
(3.22)
Procura-se minimizar Φ em C e encontrar o valor:
(3.23)
Aplicando a desigualdade de Young (ver anexo A.10) para a função indicadora I
c
(x)
(ver anexo A.6) e substituindo o vetor y por –F(x) na sua função conjugada I
c
*
(y) (ver
anexo A.4), obtém-se:
I
C
(x) + I*
C
(–F(x)) ≥ 〈 x, –F(x) 〉 para todo x em R
n
I
C
(x) + I*
C
(–F(x)) - 〈 x, –F(x) 〉 ≥ 0
I
C
(x) + I*
C
(–F(x)) + 〈 x, F(x) 〉 ≥ 0
Da definição da função indicadora, I
C
(x) ≥ 0, logo:
(3.24)
0))(()(,)(
*
≥−+=Φ xFIxFxx
C
)(inf x
Cx
Φ=
∈
α
))(()(,)(
*
xFIxFxx
C
−+=Φ
)(),(max)( yxxfxG
Cy
−∇=
∈
)(min xG
Cy∈
22
Portanto, o valor ínfimo de Φ, representado por α, é:
(3.25)
Observar que da definição da função conjugada da função indicadora I
c
:
Em particular, para z ∈C, tem-se:
Então,
Logo, pela relação em (3.19), pode-se concluir que a função Φ tem condições de se
chamar uma função gap.
Agora, considere C ⊆ R
n
um conjunto não vazio, convexo e fechado, f: R
n
→ R ∪ {∞}
uma função convexa, inferiormente semicontínua e diferenciável numa vizinhança
aberta de C. Seja g: R
n
→ R ∪ {∞} a extensão convexa da função f restrita a C
sobre R
n
definida por:
g(x) = f(x) + I
C
(x) (3.26)
Seja h: C→ R ∪ {∞} definida por:
h(x) = f(x) + g*(∇ f(x) – F(x)) + 〈 x, F(x) - ∇ f(x) 〉
(3.27)
Queremos encontrar:
(3.28)
Tomando a desigualdade de Young para g, temos:
g(x) + g*(y) ≥ 〈x,y〉 para todo x,y em R
n
(3.29)
Substituindo g(x) por seu valor em (3.26), tem-se:
f(x) + I
C
(x)+ g*(y) - 〈x,y〉 ≥ 0 (3.30)
0)(inf =Φ=
∈
x
Cx
α
)(inf xh
Cx∈
=
β
[
]
)()(,sup)(,)( zIxFzxFxx
C
Rz
n
−−+=Φ
∈
[
]
)(,sup)(,)( xFzxFxx
n
Rz
−+=Φ
∈
[
]
)(,)(,sup)( xFzxFxx
n
Rz
−=Φ
∈
[
]
))((sup)( zxxFx
n
Rz
−=Φ
∈
[
]
[
]
)(,sup)()(,sup xFzzIxFz
nn
Rz
C
Rz
−=−−
∈∈
23
Tomando y = ∇ f(x) – F(x) aplicado em (3.29)
f(x) + I
C
(x)+ g*(∇ f(x) – F(x)) - 〈 x, ∇ f(x) – F(x) 〉 ≥ 0 (3.31)
Para x ∈ C, I
C
(x) ≥ 0. Logo:
h(x) = f(x) + g*(∇ f(x) – F(x)) + 〈 x, F(x)-∇ f(x) 〉 ≥ 0 (3.32)
(3.33)
Note que quando f ≡ 0, a função h(x) torna-se igual à função Φ(x).
Após as relações desenvolvidas de (3.26)-(3.33), podem ser formuladas as
seguintes propriedades:
Propriedade 1 - Seja C ⊆ R
n
um conjunto não vazio, convexo e compacto, F: C → R
n
uma função contínua e Φ definida em (3.22). Então:
(i) Φ é não negativa, inferiormente semicontínua e limitada em C;
(ii) Φ (x) = 0 somente se x é uma solução de (3.1);
(iii) e existe um x* em C que minimiza Φ em C.
Prova: ver Auchmuty (1989).
A propriedade a seguir relaxa a hipótese de compacidade do conjunto viável C.
Propriedade 2 - Seja C ⊆ R
n
um conjunto não vazio, convexo e fechado, F:C→R
n
uma função contínua e Φ definida em (3.22). Então:
(i) Φ é inferiormente semicontínua e não negativa em C;
(ii) Φ (x) = 0 se e somente se x é solução do PDV(F, C) definido em (3.1).
Prova: ver Auchmuty (1989).
Propriedade 3 – Seja C ⊆ R
n
um conjunto não vazio, convexo e fechado; F: C → R
n
uma função contínua. Seja f : R
n
→ R ∪ {∞} uma função convexa inferiormente
0)(inf ==
∈
xh
Cx
β
0)(inf =Φ=
∈
x
Cx
α
24
semicontínua que é continuamente diferenciável numa vizinhança aberta de C e
g(x), h(x) e β definidos em 3.26, 3.27 e 3.28 respectivamente. Então:
(i) h é não negativa e inferiormente semicontínua em C;
(ii) h (x) = 0 se e somente se x é uma solução de (3.1);
(iii) Quando C é compacto, então β = 0 e h atinge seu ínfimo em C na solução
do PDV(F, C) definido em (3.1).
Prova: ver Auchmuty (1989).
3.3.2 – Formulação do Ponto Sela
Seja C ⊆ R
n
um conjunto não vazio, convexo e fechado, f: R
n
→ R ∪ {∞} uma função
convexa inferiormente semicontínua, contínua e diferenciável numa vizinhança
aberta de C e F: C → R
n
uma função contínua dada. Segundo Auchmuty (1989),
quando existe um princípio variacional baseado nos funcionais das formas (3.22) e
(3.27) há também modelos minimax e princípios variacionais duais. Desde que h
generaliza Φ, deve-se concentrar no problema (3.28).
Definimos L: C x C → (-∞, ∞) por
(3.34)
(3.35)
De (3.27), pode-se reescrever h(x) como:
h(x) = f(x) + (x, F(x) - ∇ f(x)) + [g*(∇ f(x) – F(x)) ] (3.36)
Se ∇ f(x) – F(x) = w, pode-se reescrever (3.35) e (3.36):
(3.37)
h(x) = f(x) + (x, F(x) - ∇ f(x)) + [g*(w)] (3.38)
Da definição de função conjugada para uma função estendida g, tem-se:
)()(,)()(),( xfxFyxyfxfyxL ∇−−+−=
)()(,)()(,)()(),( xfxFyxfxFxyfxfyxL ∇−−∇−+−=
[
]
)()()(,)()(,)(),( yfxFxfyxfxFxxfyxL −−∇+∇−+=
[
]
)(,)()(,)(),( yfwyxfxFxxfyxL −+∇−+=
)()(,)()(),( xFxfyxyfxfyxL −∇−−−=
25
(3.39)
(3.40)
Sabendo que I
C
(y) ≥ 0, temos que:
(3.41)
Logo (3.42)
e (3.43)
Escrevendo a função objetivo do problema dual H: C → [-∞, ∞) definida por:
(3.44)
o problema dual deverá calcular o supremo de H em C, encontrando:
(3.45)
Um ponto é chamado de ponto de sela de L em C x C se
(3.46)
para todo 〈x,y〉 em C x C.
Resumindo: O resultado dos cálculos acima diz que se a função Lagrangeana, L,
dada em (3.34) tem um ponto de sela em
)
ˆ
,
ˆ
( yx
, então
x
ˆ
será uma solução do
problema de desigualdade variacional original, PDV (F, C).
3.3.3 – Função Gap Primal e suas propriedades
Um caso especial das funções de mérito resulta de deixarmos f ≡ 0 em (3.34). A
função correspondente é:
(3.47)
É conhecido que, no contexto da teoria dos jogos não-cooperativos, esta função foi
estudada por Zuhovickii et al. (1969a, 1969b, 1970, 1973) (apud Larsson e
[
]
))()((,sup)(* yIyfwywg
C
Ry
n
+−=
∈
[
]
)(,))()((,sup)(* yfwyyIyfwywg
C
Ry
n
−≥+−=
∈
[
]
)(,sup)(* ygwywg
n
Ry
−=
∈
),(sup)( yxLxh
Cy∈
=
),(inf)( yxLyH
Cx∈
=
),(supinf yxL
Cy
Cx
∈
∈
=
β
)(sup yH
Cy∈
=
β
)
ˆ
,()
ˆ
,
ˆ
(),
ˆ
( yxLyxLyxL
≤
≤
)
ˆ
,
ˆ
( yx
)(,),( xFyxyxL −−−=
yxxFyxL
T
−= ,)(),(
26
Patriksson (1994)). Eles mostraram que, sob a hipótese de F ser monótona, a
solução para o problema de equilíbrio é encontrada buscando um ponto de sela para
L, isto é, um ponto (x
*
, y
*
) tal que:
(3.48)
Explorando a igualdade em (3.48), foi obtida a função Gap Primal G : C→R ∪ {∞},
definida por:
(3.49)
Aqui o conjunto C ⊆ R
n
é não vazio, convexo e fechado.
Esta função tem sido muito estudada em vários contextos. Na seguinte seção se
descrevem algumas de suas propriedades mais importantes.
- Propriedades:
Seja um conjunto C ⊆ R
n
não vazio, convexo e fechado, f: R
n
→ R ∪ {∞} uma função
convexa inferiormente semicontínua e diferenciável numa vizinhança aberta de C, F:
C → R
n
uma função contínua dada e Ω o conjunto solução do PDV (F, C). Para
qualquer x ∈ C, seja Y(x) o conjunto de soluções ótimas, possivelmente vazio, para
o problema de otimização definido em (3.49).
i) G é uma função Gap;
ii) G é inferiormente semicontínua em C;
iii) Se C é limitado e F ∈ C
1
em C, então G é Lipschitz contínua em C, isto é, existe
uma constante M
G
≥ 0 tal que para todo x
1
,
x
2
∈ C, tem-se:
||G(x
1
) – G(x
2
)|| ≤ M
G
|| x
1
–x
2
||;
iv) Se F ∈ C
1
em C, então G é diferenciável em x ∈ C se Y(x) = {y(x)}, com
∇G(x) = F(x) + 〈∇F(x),x-y(x)〉;
v) Se F ∈ C
1
em C, e é monótona em C, então se x ∉Ω e Y(x) = {y(x)}, a direção d =
y(x) – x é uma direção de descida viável com respeito a G em x, e a derivada
direcional G’(x;d) = 〈∇G(x),d〉 ≤ - G(x);
vi) G é convexo em C se 〈F(x), x〉 é convexo em C e cada componente de F é
côncavo em C;
yxxFxG
Cy
−=
∈
),(sup)(
)()(supinf*)*(*)()()(infsup yxxFyxxFyxxF
T
Xy
Xx
TT
Xx
Xy
−=−=−
∈
∈∈
∈
27
vii) x ∈ Ω ⇔x ∈ Y(x) é um ponto fixo característico de Ω;
viii) Sob as condições de F em (v), x ∈ Ω ⇔ G’(x, y – x) ≥ 0 ∀y ∈ C é um ponto
estacionário característico de Ω.
Prova: Em Larsson e Patriksson (1994).
3.3.4- Função Gap Dual e suas propriedades
Explorando a igualdade em (3.48), foi obtida a função Gap Dual g: C→R ∪ {- ∞},
definida por:
(3.50)
No caso geral de , (f ≠ 0), a formulação do dual
não constitui propriamente uma função gap. Para o caso de f ≡ 0, contudo, é
possível mostrar que g é uma função gap, providenciando que a função F seja
pseudomonótona em C, isto é,
F(y)
T
(x – y) ≥ 0 ⇒ F(x)
T
(x – y) ≥ 0 ∀x,y ∈ C.
- Propriedades:
Seja F : C →R
n
pseudomonótona em C. Também, para qualquer y ∈ C, seja X(y) o
conjunto (possivelmente vazio) de soluções ótimas para o problema de otimização
que se define em (3.50).
i) g é uma função gap;
ii) g é côncava em C;
iii) Se F ∈ C
1
em C, então g é diferenciável em y ∈ C se X(y) = {x(y)}, com
∇g(y) = - ∇F(x(y)).
Prova: Em Larsson e Patriksson (1994).
yxxFyg
Cx
−=
∈
),(inf)(
)()(,)()(),( xFxfyxyfxfyxL −∇−−−=
28
3.3.5 – Exemplo de uma função Gap para um programa Convexo
Considere o seguinte problema, estudado em Hearn (1981), de programação
convexa:
min f(x) s. a. x ∈ S = { x / Ax = b, x ≥ 0}, (3.51)
para f convexo e continuamente diferenciável e S um
polítopo não-vazio e compacto,
tal que existe um x* que resolva (3.51).
A proposta desenvolvida neste trabalho é resolver (3.51) indiretamente,
considerando uma função auxiliar gap, G(x), definida por:
G(x)
(3.52)
Dadas as seguintes propriedades da função G(x):
G(x) ≥ 0 para todo x ∈ S, (3.53)
G(x*) = 0 se e somente se x* resolve (3.51) (3.54)
Pretende-se mostrar no decorrer desta seção que a solução do problema (3.51)
pode ser obtida através do problema auxiliar ( que tem valor ótimo igual a zero):
min G(x), s. a. x ∈S (3.55)
A função Lagrangiana do problema (3.51) é dada por L:Ω x R
l
x R
m
→ R
L (x, λ, µ) = f(x) +
L (x, λ, µ) = f(x) + , que aplicada ao problema (3.51) é:
h(x) = b – Ax = 0
g(x) = - x ≤ 0
λ ∈ R
µ ∈ R
+
L (x, λ, µ) = f(x) +
L (x, λ, µ) = f(x) - 〈 λ, Ax – b 〉 - 〈 µ, x 〉
)(),(max yxxf
Sy
−∇=
∈
)(),(min xyxf
Sy
−∇−=
∈
∑∑
==
+
m
j
jj
l
i
ii
xgxh
11
)()(
µλ
)(,)(, xgxh
µλ
+
xAxb −+− ,,
µλ
29
Seja D
Quando x ∈ D, para todo (λ, µ) ∈ R
l
x R
m
+
⇒ L (x, λ, µ) ≤ f(x)
Então, o Problema (3.51) pode ser reescrito como:
sujeito a x ∈ D. (3.56)
O dual do problema (3.56) é:
(3.57)
Onde ∆ = {(λ,µ) ∈ R
l
x R
m
+
/ inf L (x, λ, µ) > - ∞ }
Como é conhecido que o , quando a função objetivo e as funções de
igualdade e desigualdade são diferenciáveis, ocorre quando ∇
x
L (x, λ, µ) = 0, logo
podemos reescrever o problema dual como:
(3.58)
sujeito a ∇
x
L (x, λ, µ) = 0
µ≥ 0
onde L (x, λ, µ) = f(x) - 〈 λ,(Ax-b) 〉 - 〈 µ,x 〉, logo
∇
x
L (x, λ, µ) = ∇ f(x) - 〈 A,λ 〉 - µ
Definimos
D(x) = {(λ, µ)/∇
x
L (x, λ, µ) = 0, µ ≥ 0} ( conjunto viável dual)
( valor ótimo da função Lagrangiana para cada x em D(x))
Estudando a função Gap de um programa convexo, Hearn (1981) conseguiu
destacar duas importantes propriedades:
Propriedade 1:
x ∈ S implica D(x) ≠ ∅
),,(max
),,(
µλ
µλ
xL
x
),,(max)(
)(),(
µλ
µλ
xLxd
xD∈
=
+
∈
),,(supmin
),(
µλ
µλ
xL
ml
xRR
+∞<Ω∈=
+
∈
),,(sup
),(
µλ
µλ
xLx
ml
xRR
{
}
),,(inf max
,
µ
λ
µλ
xL
x Ω∈
∆∈
),,(inf
µλ
xL
x Ω∈
30
Prova em Hearn (1981).
Propriedade 2:
d (x) = f (x) – G (x) ∀ x∈S.
Prova em Hearn (1981).
Logo, se for encontrado o valor G(x) = 0, tem-se d(x) = f(x), o que faz concluir que a
função f(x) já chegou ao seu valor ótimo, portanto se resolve o problema.
3.4 – Método de Penalização
Nesta seção são apresentadas opções da estratégia de taxar para inibir ocorrências
pouco desejadas, isto é, aplicar uma penalidade. Este artifício em programação
matemática tem por objetivo obter uma configuração mais simples do modelo em
estudo, conseguindo, assim, calcular soluções com os vários métodos conhecidos
que existem na literatura. Como se viu no capítulo 2, uma outra opção de aplicar
penalidade foi a de colocar em evidência os casos de desequilíbrio na entrega de
gás natural, que tem a dificuldade de armazenagem, portanto, qualquer excesso ou
falta do gás ocasiona maiores custos.
3.4.1 – O conceito de Função Penalizada
Segundo Bazaraa et al. (1993), as funções penalizadas surgem como uma
alternativa de transformar um problema com diversas restrições em um outro
problema ou família de problemas sem restrições. A estratégia usada é colocar as
restrições dentro da função objetivo através de um parâmetro de penalização que
inibirá qualquer violação das restrições. Para motivar o estudo, considere o seguinte
problema de minimização com restrição de igualdade:
minimizar f(x)
sujeito a: h(x) = 0 (3.59)
Seja µ > 0 um número bem grande. O problema (6.1) pode ser reescrito como:
minimizar f(x) + µ h
2
(x)
31
sujeito a: x ∈ R
n
(3.60)
Observa-se que a solução ótima do problema (3.60) deve ter h
2
(x) perto de zero,
pois se isso não se realiza, o fator da penalidade µ imposta aumentará o valor da
função, o que não é desejado para obter o mínimo do problema original.
Agora considere o seguinte problema com restrição de desigualdade:
minimizar f(x)
sujeito a: g(x) ≤ 0 (3.61)
A forma f(x) + µ g
2
(x) não é apropriada, pois esta penalidade pode deixar a restrição
assumir valores g(x) > 0, isto é, trabalhar com pontos que não são viáveis. Neste
caso, o artifício para a penalidade deve apenas taxar o fato de x não ser viável.
Assim, o problema (3.61) pode ser reescrito como:
minimizar f(x) + µ Máximo{0, g(x)}
sujeito a: x ∈ R
n
(3.62)
Se g(x) ≤ 0, então o Máximo {0, g(x)} = 0 e a parcela da penalidade assumirá valor
zero. Se g(x) > 0 então o Máximo {0, g(x)} > 0 e a penalidade µ g(x) aparecerá com
força no problema penalizado.
Observe que no ponto x onde g(x) = 0, a função objetivo penalizada pode não ser
diferenciável embora f e g sejam diferenciáveis. Como a diferenciabilidade é
desejável, a parcela da penalidade em (3.62) será formulada por:
µ [Máximo {0, g(x)}]
2
.
Em geral, uma função penalizada deve colocar uma punição positiva alta, para os
pontos que não viáveis, e nenhuma punição para os pontos viáveis. Se as restrições
são da forma g
i
(x) ≤ 0 para i = 1,...,m e h
j
(x) = 0 para j = 1,..., l, então uma função
penalizada α(x), Bazaraa et al. (1993) definem por:
(3.63)
[ ]
[ ]
)()()(
11
xhxgx
j
l
j
i
m
i
∑∑
==
+=
ψφα
32
Onde φ e ψ são funções contínuas satisfazendo o seguinte:
φ(y) = 0 se y ≤ 0 e φ(y) > 0 se y > 0
ψ(y) = 0 se y = 0 e ψ(y) > 0 se y ≠ 0 (3.64)
Tipicamente, φ e ψ são da forma:
φ(y) = [Máximo{0, y}]
p
ψ(y) = |y|
p
Onde p é um inteiro positivo. Assim, a função penalizada α(x) é usualmente da
forma:
(3.65)
3.4.2 - Penalidade Exata e Inexata para modelos de Otimização com
Restrições
Um método analítico e algorítmico que resolve um modelo de otimização não linear
com restrições
minimizar f(x)
s.a. x
∈
X ,
(3.66)
onde
{
}
kkxgmixhxX
ji
n
,...2,1 ,0)( ;,...,2,1 ,0)(| =≤==ℜ∈= é a penalização exata,
Bazaraa et al. (1993). A idéia do método consiste em resolver um modelo ou uma
seqüência de modelos de minimização sem restrições de igualdade, nem
desigualdade, onde à função de custo é adicionado um termo que penaliza
fortemente as diferenças que possam existir com as restrições do modelo original.
Um exemplo importante de função de custo (função objetivo) penalizada ou função
auxiliar quadrática é:
[ ]
[ ]
++=
∑∑
=
+
=
k
j
j
m
i
xgxhxfxf
1
2
1
2
)()()()(
α
α
(3.67)
sendo
{
}
.)(,0max)( xgxg
jj
=
+
Observa-se que associado ao termo da penalização
está o parâmetro α > 0, que determinará a severidade da multa pela violação das
{ }
[ ]
∑∑
==
+=
l
j
p
j
p
m
i
i
xhxgmáximox
11
)()(,0 )(
α
33
restrições. Observe que a função penalizada herda a não diferenciabilidade
)(xg
j
+
,
porque é definida como o valor máximo entre dois valores.
Agora, considere a notação
α
x
para o vetor solução e
θ
(
α
) como o valor ótimo do
modelo auxiliar:
[
]
[
]
{
}
∑∑
=
+
=
∈
++=
k
j
j
m
i
i
Rx
xgxhxfxf
n
1
2
1
2
)()( )()( min
α
α
=:
θ
(
α
), (3.68)
Deve-se notar que as soluções de
θ
(
α
) de cada problema penalizado, possivelmente
sejam pontos não viáveis do modelo de otimização original (3.66) para valores finitos
de
α
, mas ao aumentar o valor de
α
ficaria muito mais garantido que as restrições do
modelo original fossem satisfeitas; configurando-se assim, em uma estratégia de
penalização externa, como é conhecida na literatura.
A seguir apresenta-se o resultado mais relevante sobre a abordagem de problemas
de otimização com restrições via modelo penalizado.
Em Bazaraa et al. (1993) o teorema 9.2.2 , pagina 367, garante que:
(i)
},...,2 ,1 ,0)( ;,...,2 ,1 ,0)(||)(inf{ kkxgmixhRxxf
ji
n
=≤==∈
=
)(lim
α
θ
α
∞→
;
(ii)
o
xx =
∞→
α
α
lim
de qualquer seqüência convergente de soluções
{
}
α
x
do modelo
penalizado é solução do modelo original;
(iii) o termo
[
]
[
]
{
}
0)()(
1
2
1
2
→+
∑∑
=
+
=
k
j
j
m
i
i
xgxh
α
quando
∞
→
α
sob as condições de continuidade das funções f, h
i
e g
j
para i=1, 2,...,m; j=1, 2, ...k,
da existência de pontos viáveis para o modelo na relação (3.66), continuidade de f
α
,
a função auxiliar, e a existência de solução do modelo sem restrições na relação
(3.68).
Porém, a minimização da função de custo penalizada, ou da função auxiliar
quadrática, apresentada na relação (3.68) não é exata. Uma função de custo
penalizada é dita exata se existe certo α
0
> 0 tal que uma solução da minimização
dessa função de custo penalizada é solução do problema original (3.66) para
quaisquer valores de α ≥ α
0
. Na década de 70, foi estudada uma versão diferente de
função de custo penalizada:
34
[
]
{
}
∑
∑
=
+
=
++=
k
j
j
m
i
i
xgxhxfxF
11
)()()()(
α
α
(3.69)
Usando esta função, mostrou-se que para certas funções, uma solução local x* de
(3.66) é solução local de:
(3.70)
sob condições intimamente ligadas aos multiplicadores de Lagrange das restrições
de (3.66), proporcionando um certo valor de
α
a partir do qual x* é solução local do
modelo (3.70) para todo α ≥
α
. A seguir se transcreve o teorema 9.3.1 de Bazaraa
et al. (1993), página 374, que estabelece a grandeza da penalidade
α
a partir da
qual o modelo penalizado tem a mesma solução local que o modelo original, sob
condições de convexidade para as desigualdades e linearidade para as restrições de
igualdade.
Teorema 1 - Considere o modelo
minimizar f(x)
s.a.
;,...,2,1 ,0)(
mixh
i
==
kkxg
j
,...2,1 ,0)( =≤ .
Sejam as funções g
j
convexas e as h
i
afins e diferenciáveis. Supor que o vetor (x°, v,
u), onde x° é uma solução viável para o modelo, e os v
i
, u
j
são os multiplicadores
de
Lagrange associados às restrições de igualdade e desigualdade respectivamente,
satisfazendo as condições de Karush-Kuhn-Tucker:
(3.71)
Então, para
α
= máximo{ u
j ,
|v
i
|
}, o ponto x° é solução local do modelo (3.70) para
todo
αα
≥
, isto é, minimiza a função de custo penalizado
α
F , de modo exato, para
cada
αα
≥
.
3.4.3 - Penalidade Exata e Inexata para problemas de otimização com
uma Desigualdade Variacional
0 321
0)(
0)()()(
11
≥=
=°
=°∇+°∇+°∇
∑∑
==
j
jjj
k
j
jj
m
i
ii
u,...k , , para j
xguativogcadapara
xguxhvxf
[
]
{
}
∑
∑
=
+
=
∈
++=
k
j
j
m
i
i
Rx
xgxhxfxFimizar
n
11
)()()()(min
α
α
35
Um problema de programação matemática onde uma restrição é representada por
uma desigualdade variacional, é conhecido na literatura como um Problema
Generalizado de Programação Matemática de Dois Níveis (GPMDN) (Marcotte e Zhu
(1996) e Morales (1997)) ou também problema de programação matemática com
uma restrição de equilíbrio. Assim, um GPMDN é escrito como:
Min f(x, y)
s.a. y solução de
(3.72)
〈
F(x,y) (z – y)
〉
≥ 0
∀
z
∈
C(x),
Para este problema, se assume ter um conjunto S de definição convexo, F: S →R
n
continuamente diferenciável e monótono sobre C(x) para cada x viável.
Nesta seção, se apresenta primeiro uma nova versão que define o problema com
uma desigualdade variacional, como restrição, e termina com a proposta de cálculo
de solução de um problema GPMDN mediante o algoritmo que penaliza a função
gap associada ao PDV(F(x, ⋅), C(x)).
Inicia-se esta seção com uma alternativa para a compreensão de uma desigualdade
variacional e sua estreita relação com a condição de otimalidade de primeira ordem
do modelo de otimização não linear (PNL) em (3.66), onde se relaxa a condição de
diferenciabilidade das funções. De fato, lembrar o modelo:
minimizar f(x)
s.a. x
∈
X ,
onde X será trabalhado sem especificação detalhada, ou seja, X será um conjunto
abstrato. Este problema, teoricamente, pode-se resolver trabalhando com um
problema de minimização irrestrito, onde a função objetivo é f + I
X
, com I
X
sendo a
função indicadora do conjunto de restrições X. Caso o conjunto X seja um conjunto
convexo, então a funçãondicadora, I
X
, é uma função convexa.
Nessas condições se conhece o seguinte resultado.
Teorema 2 - Seja X um conjunto convexo e fechado e f uma função convexa, então
as seguintes afirmações são equivalentes quando
X
x
∈
:
(i)
x
minimiza f sobre o conjunto X;
(ii)
)]([0 xIf
X
+∂∈ , isto é, o vetor nulo pertence ao operador subdiferencial (ver
anexo A.8) da função f + I
X.
36
Prova: Teorema VII 1.1.1, p. 293 em Hiriart-Urruty e Lemachal (1993).
Assumindo que a função f é convexa e diferenciável, obtém-se:
)()()()()]([ xNxfxIxfxIf
XXX
+∇=∂+∇=+∂
onde
)(
xN
X
é o cone normal ao conjunto X no ponto
x
. Usando a definição de cone
normal (ver anexo A.12), pode-se reescrever a condição (ii) do Teorema 2
como
)()(0 xNxf
X
+∇∈ . Portanto, resolver o problema de PNL consiste em
encontrar um
X
x
∈
tal que )()( xNxf
X
∈∇− ,
0),( ≥〉−〈∇ xyxf para todo
Xy
∈
(3.73)
Assim escrito não é difícil reconhecer as condições geométricas de otimalidade de
primeira ordem de um PNL. O problema de desigualdade variacional é uma
generalização de (3.72), onde se substitui
)(
xf
∇
por uma função F: R
n
→R
n
. Assim
definido,
achar
X
x
∈
de modo que 0),( ≥〉−〈 xyxF para todo
Xy
∈
. (3.74)
Para resolver problemas generalizados de programação matemática, onde uma das
suas restrições é descrita por uma desigualdade variacional, um PDV (F, C) definido
em (3.73), Marcotte e Zhu (1996) trabalharam usando as propriedades de uma
função regularizada gap associada a PDV(F,C) para construir métodos de descida,
definida por ),(max)( zxxG
Xz
φ
α
∈
= sendo
2
||||)
2
(),(),( zxzxxFzx −−〉−〈=
α
φ
, α positivo (3.75)
Marcotte e Dussault (1989) utilizaram a função gap G
α
quando α=0, a mesma
apresentada na seção 3.3 resultando uma função gap linear para um problema
generalizado de programação linear. Fukushima (1992) e Trujillo (2003) trabalharam
e mostraram que ),(max)( zxxG
Xz
φ
α
∈
= é uma função gap quadrática diferenciável
para α > 0 . A função
φ
(x, z) é côncava na variável z, e pelas propriedades de uma
função gap, ver seção 3.3 e 3.4, G
α
é não negativa no domínio do PDV(F,C) e é
igual a 0 quando se está na solução do mesmo.
Portanto, no problema GPMDN em (3.72), o segundo nível do modelo, a
desigualdade variacional denotada por PDV(F(x,
⋅
), C(x)), possui uma função gap
parametrizada sendo a variável x o seu parâmetro:
37
2
)()(
||||)
2
(),,(max),,(max),( zyzyyxFzyxyxG
xYzxYz
−−〉−〈==
∈∈
α
φ
α
Em conseqüência pode-se escrever o problema GPMDN (3.72) como segue:
Min
x
∈
X, y
∈
Y(x)
f(x, y)
s.a. (3.76)
G
α
(x,y)= 0.
Para calcular a solução do problema em (3.76), Marcotte e Zhu (1996) propuseram a
aproximação mediante a família de problemas penalizados a seguir:
),(),(),,( min
,
yxGyxfyxQ
Syx
αα
µµ
+=
∈
(3.77)
com µ um número positivo.
Marcotte e Zhu (1996) mostram que sob as hipóteses iniciais do problema, (3.77) é
um problema de programação não linear com restrições convexas e função objetivo
diferenciável para α positivo. Adicionalmente provaram que a função penalizada
formulada em (3.77) é uma função de penalidade exata sob as seguintes condições:
domínio S poliédrico, a função F fortemente monótona com respeito à variável y e
com seu jacobiano,
F
y
∇
, Lipschitz contínuo sobre S. Em geral é uma função de
penalidade não exata.
38
Capítulo 4
PRÉ-PROCESSAMENTO
Como foi visto no capítulo 2, o modelo apresentado é um problema de programação
matemática em dois níveis misto. O conjunto de soluções ótimas, associado com um
problema de programação em dois níveis é determinado implicitamente por uma
série de dois problemas de otimização, que devem ser resolvidos em uma seqüência
predeterminada: Temos um líder que seleciona sua decisão primeiramente e o
seguidor que responde a esta decisão. Entretanto, quando a função objetivo e/ou
algumas restrições do problema não são lineares (o que acontece no modelo do
capítulo 2 em estudo), resolver uma série de dois problemas de otimização, neste
caso particular minimização, torna-se uma tarefa complexa.
A proposta é utilizar técnicas algébricas conhecidas para manipular o modelo do
problema, tornando-o mais apropriado para os métodos de solução existentes para
modelos de programação matemática de apenas um nível, respeitando as condições
do problema original, de modo a manter a equivalência em relação às soluções, o
que na literatura é chamado de pré-processamento (Morales e Ríos-Mercado, 2002).
Sendo assim, o passo inicial é transformar o modelo matemático do problema em
outro equivalente, de forma que tenha melhores propriedades computacionais para
ser resolvido numericamente por métodos conhecidos.
4.1 – Reformulação do Modelo do Problema de Dois Níveis
Kalashnikov e Ríos-Mercado (2006) apresentam uma proposta de reformulação do
modelo mostrado no capítulo 2. A proposta de transformar o problema de dois níveis
em um problema de apenas um nível será colocada a seguir:
Em primeiro lugar, cabe observar que a função objetivo do segundo nível, |z|, não é
diferenciável. Isto impede, por exemplo, o cálculo do seu gradiente e a aplicação de
técnicas que necessitem dele. Para evitar este tipo de dificuldade, resolveu-se
regularizar esta função transformando-a em z
2
, que é uma função diferenciável.
39
Deve-se atentar, também, para o fato de que, quando os desequilíbrios do último dia
em todas as zonas são não-positivos ou positivos e grandes o suficiente, o conjunto
solução do problema de dois níveis (2.1a) – (2.2l) não mudará se a restrição ligada à
variável inteira q em (2.2h) for movida do segundo para o primeiro nível (Dempe et
al. (2005)). Sendo q a única variável inteira, basta fixá-la no modelo original e rodá-
lo para q=0 e para q=1, comparando as respostas e assumindo como ótima a de
menor valor. Assim, o modelo é definido por:
minimizar h
1
(x,s,y,u,v,z) = z
(x,s)
s.a.
(2.1b)-(2.1e)
minimizar h
3
(x,s,y,u,v,z) = z
2
(y,u,v)
s.a.
(2.2b)-(2.2k)
q = 0
minimizar h
1
(x,s,y,u,v,z) = z
(x,s)
s.a.
(2.1b)-(2.1e)
minimizar h
3
(x,s,y,u,v,z) = z
2
(y,u,v)
s.a.
(2.2b)-(2.2k)
q = 1
Estes são sistemas hierárquicos em dois níveis no espaço Euclidiano, onde a
tomada de decisão do primeiro nível controla um vetor de variáveis w
1
= (x,s) ∈ R
NP
x
R
NP
, e o segundo nível controla um vetor de variáveis w
2
=(y,u,v) ∈ R
P
x R
P(P-1)/2
x
R
P(P-1)/2
, e uma variável binária q = β. O líder toma sua decisão primeiro,
considerando que ela influenciará a decisão do seguidor. Definindo a função
objetivo como a função:
∑
∑
∑
<∈ <
+
−−+−=
jiji
ijijij
Ji jiji
ijijiiii
uefvbyyrvuyF
:),(:),(
2
)1(])([),,(
δ
40
E o conjunto das soluções do seguidor:
W
2
β
= W
2
β
(x,s) = {w
2
β
= (y,u,v) : (2.2b)-(2.2h), (2.2j)-(2.2k), fazendo q = β}
Então, w
2
β
= w
2
β
(x,s) = (y,u,v) é a reação ou resposta do seguidor para a decisão
do líder w
1
= (x,s), e também uma solução do problema de equilíbrio representado
pela seguinte desigualdade variacional (na seção 3.4.3 a condição de otimalidade de
primeira ordem do problema do seguidor):
(4.1)
Notar que esta função F(y,u,v) é o negativo da função objetivo original, onde se
incluiu um termo com a penalidade contratual do emissor. Desta forma se pode
obter os seguintes problemas generalizados de programação matemática em dois
níveis equivalentes, ou GPMDN(β), para β = 0,1, na forma de maximização:
s.a. w
2
β
solução
Onde o líder está implicitamente restrito ao conjunto W
1
de vetores (x,s) tais, que o
conjunto de restrições W
2
(x,s)
β
do segundo nível é não vazio. A solução ótima deste
problema pode ser obtida da resolução de ambos os problemas GPMDN(0) e
GPMDN(1).
Agora, a desigualdade variacional do segundo nível será reformulada como uma
função gap associada para resolver o problema de segundo nível. Usamos então
esta função “gap” como um termo de penalidade para o problema de primeiro nível.
Desde que a função gap caracterize que o segundo nível é não-negativo sobre o
domínio viável, o termo de penalidade se reduz a uma forma muito simples. Esta é a
função gap associada com a desigualdade variacional (4.1) do segundo nível:
(4.2)
Onde
),,,(max),(
221
),(
21
2
2
wwwwwG
sxWw
βββ
α
φ
β
∈
=
,
2
1
),()(),,(
2
222222221
βββββ
αφ
wwwwwFwFwww −−−∇−=
T
jijiijjijiijijJiiii
befyrFcom })(,))1((,])(2{[(
:),(:),(
<<∈+
−−−=∇
δ
),,(,0),()(
222222
sxWwtodoparawwwFwF
ββββ
∈≥−∇
),,,(max:)(
),(),,,(),(
,*
2
1
vuyFFGPMDN
sxWvuyWsx
β
β
β
∈∈
=
),,(,0),()(
222222
sxWwtodoparawwwFwF
ββββ
∈≥−∇
41
e α é um número não-negativo.
A função
φ
é côncava em w
2
(fortemente côncava se α é positivo). Também, G
α
β
é
não-negativo sobre:
Se e somente se
β
2
w é uma solução para a desigualdade variacional do segundo
nível parametrizada por (x,s), a desigualdade variacional (4.1) pode ser reescrita
com a seguinte equação não-linear:
onde ),(
21
β
α
wwp é qualquer solução de (4.2). Finalmente isto nos conduz à
reformulação do GPMDN(β) como um problema de programação matemática
padrão:
(4.3)
sujeito a
=),(
21
ββ
α
wwG
0
q =
β
Para implementar o problema PR1(β), foi feita uma reformulação do problema
segundo Kalashnikov e Ríos-Mercado (2006):
(4.4)
Onde µ é um número positivo fixo e o subconjunto C
β
é dado por:
C
β
= {( w
1
,w
2
β
) : (2.1b) – (2.1e), (2.2b) – (2.2h), (2.2j) – (2.2k), q = β}.
O Lema seguinte garantirá a existência de solução da família de problemas
penalizados ),,(
21
µ
β
α
wwQ :
Lema 1: O subconjunto C
β
é compacto para cada β = 0,1.
Prova: Em Kalashnikov e Ríos-Mercado (2006).
0),(),(
21121
11
=×∪=
∈
ββ
α
β
β
wwGewWwS
Ww
,0)),(,,(),,(max),(
2121221
),(
21
2
2
===
∈
β
α
ββββ
α
φφ
β
wwpwwwwwwwG
sxWw
),,,(max:)(1
),(),,(,),(
22
11
vuyFPR
sxWvuywWsxw
ββ
β
∈=∈=
),()(min),,(min:)(2
212
),(
21
),(
2
1
2
1
ββ
α
ββ
α
µµβ
β
β
β
β
wwGwFwwQPR
CwwCww
+−=
∈∈
42
Então, o problema PR2(β) é não-linear e sua função objetivo é continuamente
diferenciável quando α é positivo. Para cada valor de µ se denota por (w
1
(µ),w
2
β
(µ))
uma solução ótima global de PR2(β) que sempre existe, devido ao fato do conjunto
viável C
β
ser compacto e a função objetivo Q
α
ser contínua. Então, um algoritmo de
penalização inexata é obtido especificando uma seqüência de valores {µ
k
},
crescentes e positivos, obtendo uma seqüência de soluções {(w
1
(µ),w
2
β
(µ))}
associadas, para cada valor de µ, segundo Marcotte e Zhu (1996).
Teorema 1: Seja {(w
1
k
, w
2
k
)}
∞
k=1
uma seqüência de iterações geradas por um
algoritmo de penalidades baseado na função Q
α
do programa PR2(β).
Então, todo ponto limite desta seqüência é uma solução do programa de dois níveis
PR1(β).
Prova: Em Kalashnikov e Ríos-Mercado (2006).
Observação: Depois de obter as duas soluções ótimas (ŵ
1
β
,ŵ
2
β
) para os problemas
PR1(β), β = 0,1, podemos achar a solução ótima do problema inicial (2.1a) – (2.1e),
(2.2a) – (2.2l) como segue:
4.2 – Reformulação das restrições do problema de dois níveis
A reformulação do modelo do problema apresentada na seção anterior foi usada
para implementar um algoritmo direto, usando o algoritmo simplex Nelder-Mead, que
busca a solução ótima do modelo do problema. Para maiores detalhes veja
Kalashnikov e Ríos-Mercado (2006).
A proposta deste trabalho é implementar um algoritmo para buscar a solução do
modelo do problema, usando programas encontrados no ambiente MATLAB,
adequando-os para o caso presente, com condições de não diferenciabilidade.
≥
=
.),,(
),()(),,(
),(
1
2
1
1
1
2
0
2
0
2
0
1
*
2
*
1
contráriocasoww
wFwFseww
ww
43
Para implementar o problema (4.4), apresentado por Kalashnikov e Ríos-Mercado
(2006), foi necessário adequar a escrita de algumas restrições para que fosse
possível utilizar a sub-rotina para otimização encontrada no MATLAB, chamada
fmincon, que encontra o mínimo de uma função de várias variáveis não linear,
sujeita a restrições de igualdade e desigualdade, onde se assume as funções duas
vezes diferenciáveis. Esta sub-rotina implementa o método seqüencial quadrático,
fazendo estimativas da hessiana do Lagrangeano do problema, e usando o método
BFGS.
Para adequar o modelo do problema (4.4) sujeito a C
β
= {( w
1
,w
2
β
) : (1b) – (1e), (2b)
– (2h), (2j) – (2k), q = β}, reescrevemos as restrições de C
β
, que não se adequavam
as entradas da sub-rotina fmincon.
As restrições modificadas (2.2e) e (2.2f) dispostas a seguir tem como dificuldade a
não diferenciabilidade, porém testes numéricos garantiram bom comportamento:
uij = min ( |xni|, max(0, xni) * max(0, -xnj)) (2.2e)
vij = min (|xnj|, max(0, -xni) * max(0, xnj)) (2.2f)
Este pré-processamento faz com que as restrições do modelo possam ser usadas
na minimização usando o fmincon.
O próximo passo é começar as experiências no ambiente MATLAB, que será tratada
no próximo capítulo.
44
Capítulo 5
EXPERIÊNCIAS NUMÉRICAS
Para avaliar e verificar o modelo do problema do transporte de gás apresentado em
(4.4), Kalashnikov e Ríos-Mercado implementaram um algoritmo de busca direta, no
qual é usado o algoritmo simplex de Nelder- Mead para resolver o problema do
primeiro nível. Para maiores detalhes veja Kalashnikov e Ríos-Mercado (2006).
A proposta deste trabalho é de implementar um novo algoritmo que utiliza os
métodos seqüencial quadrático, quase Newton com busca linear, que formam parte
da biblioteca de otimização do MATLAB, e comparar os resultados com os do artigo
estudado Kalashnikov e Ríos-Mercado, (2006).
Entretanto, antes de implementar o algoritmo para o problema do gás natural
apresentado no capítulo 2, foi conveniente trabalhar com problemas testes de
otimização (neste caso, minimização) cujos resultados numéricos sejam simples de
serem conferidos, para então trabalhar com na solução do modelo propriamente dito.
Com estes problemas testes é possível verificar os resultados numéricos,
teoricamente estabelecidos no capítulo 3. Assim, o processo de analisar e aplicar o
algoritmo desenvolvido ao problema torna-se mais confiável.
A seguir, apresentam-se os resultados de experiências numéricas realizadas para
este trabalho.
5.1 – Primeiro Problema teste
Para começar, foi escolhido um problema de apenas um nível. Este problema foi
formulado devido à facilidade de resolução.
Min f(x)
(5.1)
s. a. x ∈ C
Onde f(x) = x
2
e C = [1,3].
45
A nossa proposta é mostrar que a solução deste problema de minimização é igual à
solução do problema:
Min G(x) (5.2)
s. a. x ∈ C
Onde G(x) é uma função gap associada à desigualdade variacional, que não é mais
do que a condição necessária de primeira ordem do problema (5.1).
O primeiro passo é montar um PDV (F, C).
Se f(x) = x
2
, então ∇f(x) = 2x. Logo o PDV (F, C) pode ser montado da seguinte
forma:
Encontrar um x ∈ C que satisfaz 〈 2x, (y - x) 〉 ≥ 0, ∀ y ∈ C. (5.3)
Para resolver o PDV (F, C), precisamos encontrar x ∈C tal que:
2x (y - x) ≥ 0, para todo y ∈ [1,3],
2xy - 2x
2
≥ 0
2xy ≥ 2x
2
Observação: Apesar de não ser necessário resolver esta desigualdade, não é difícil
concluir que a solução deste PDV (F, C) é x = 1.
O segundo passo é formular a função gap associada a este problema de
desigualdade variacional que pode ser definida por:
(5.4)
Como foi visto na seção 3.2, a solução do PDV (F, C) ocorre quando G(x) = 0; e na
seção 3.3, se apresentou o comportamento da função gap G(x) ≥ 0, então, para
achar a solução do PDV (F, C) bastará resolver:
(5.5)
A solução deste problema ocorre no ponto x=1 e y=1.
)(,2max)( yxxxG
Cy
−=
∈
xyxxG
Cy
22max)(
2
−=
∈
−
∈
∈
xyx
Cy
Cx
22maxmin
2
46
Como este exemplo é pequeno, foi possível achar a sua solução sem fazer uso de
algoritmos. Desta forma, ao se implementar a resolução com um algoritmo é
possível saber se a solução está sendo calculada de maneira correta. Para isso foi
escolhido o programa MATLAB 7.0. Neste programa encontramos uma sub-rotina de
minimização pronta para uso, chamada fmincon. Ela foi preparada para ser usada
em minimizações de funções com restrições. Para usá-la, adaptamos a fórmula (5.5)
em outra equivalente:
(5.6)
Observe que se trabalha com duas minimizações.
Na primeira minimização a ser resolvida, a que está entre colchetes, a variável é y.
A sub-rotina de miminização usada, fmincon, não aceita a variável x sem valor
numérico, pois a rotina não trabalha com duas variáveis a menos que sejam
componentes de um vetor (incógnita). A alternativa usada foi aplicar a rotina
somente a função função gap em pontos discretos x ∈[1,3] e verificar a resposta do
argmin = y*.
Ao fazer este procedimento na minimização em questão, foi observado que para
diferentes valores de x o valor de y não varia, é sempre y=1, o que se esperava da
análise algébrica desta minimização, para todo valor de x∈C. Sendo assim, assume-
se o valor y* =1 e resolve-se a segunda minimização, que agora pode ser escrita
como:
(5.7)
Na minimização mais externa, a variável é x. Assim, ao final da minimização se tem
dois valores fixos x* e y*, que são os argumentos das minimizações, que
representam o valor ótimo para o problema (5.6).
Observe que a solução x*=1 e y*=1 ao ser substituída na fórmula (5.4) encontra-se o
valor gap = 0. De acordo com o que foi visto no capítulo das fundamentações
teóricas, quando o gap=0 a solução encontrada para (5.6) também é solução de
(5.3) e conseqüentemente é solução de (5.1). Isto é, o problema min-max (5.6) é
+−−
∈∈
xyx
CyCx
22minmin
2
xx
Cx
22 min
2
−
∈
47
resolvido aplicando a rotina fmincon para achar o valor da função gap em valores
discretos da variável x, dos quais se concluiu a solução final.
Seguindo o mesmo raciocínio, um segundo exemplo foi escolhido para aplicar as
técnicas utilizadas neste primeiro.
5.2 – Segundo Problema teste
Este problema foi adaptado da dissertação de Trujillo (2003).
Seja f: R
n
→R
n
uma função diferenciável cujo gradiente ∇f: R
n
→R
n
é chamado F(x)
e definido da forma: F(x) = Mx+q, onde
[M]
ii
= 4(i-1) + 1, i = 1,...,n,
[M]
ij
= [M]
ii
+ 1, i = 1,..., n-1, j = i+1,...,n,
[M]
ij
= [M]
jj
+ 1, j = 1,..., n, i = j+1,...,n,
e o vetor q = (-1,...,-1)
T
.
Seja o conjunto C dado por: C={x∈ R
n
: 0≤ x
i
≤ 50, i = 1,..., n}.
O segundo problema teste é definido por:
Min f(x)
(5.8)
s. a. x ∈ C
O primeiro passo para encontrar a solução deste problema é montar um PDV (F, C),
que pode ser escrito da seguinte forma:
Encontrar um x ∈ C que satisfaz 〈 F(x), (y - x) 〉 ≥ 0, ∀ y ∈ C
(5.9)
A solução desta desigualdade variacional encontra-se em Trujillo (2003), escrito da
forma x* = (1,0,...,0)
T
.
A proposta é mostrar que a solução do problema (5.8), assim como no primeiro
exemplo,é equivalente à solução do problema:
48
Min G(x) (5.10)
s. a. x ∈ C
Onde G(x) é uma função gap associada ao PDV (F, C) em (5.9) que pode ser
definida por:
(5.11)
Como concluímos no primeiro exemplo, para achar o resultado do PDV (F, C) basta
fazer:
(5.12)
Assim, a solução de (5.12) também é solução de (5.8).
Nossos testes foram feitos para n=2, n=4 e n=16. Mostraremos os passos de
maneira mais detalhada assumindo n=2.
Seja n = 2. A função F: R
2
→R
2
é definida por:
F(x) =
Logo (5.12) pode ser reescrita como:
Observação: Se for colocado o valor de F(x) em (5.9), tem-se:
Encontrar um x ∈ C que satisfaz:
Então, reescrevendo, tem-se:
Conhecendo o conjunto C, sabe-se que o lado direito da desigualdade é sempre
positivo. Já o lado esquerdo só será positivo se:
)(),(max)( yxxFxG
Cy
−=
∈
+−−−+++−−−+−−
∈∈
)525222(minmin
222212
2
22111211121
2
1
yyxyxxxxxyyxyxxxxx
CyCx
−−
∈∈
)(),(minmin xyxF
CyCx
Cyyyxyxxxxxyyxyxxxxx ∈∀≥−+++−−−+++−− ,0525222
222212
2
22111211121
2
1
2
2
221121
2
12222111211
522522 xxxxxxxxyyxyxyyxyx −++−+≥−++−+
2
22121
2
1212211
54)152()12( xxxxxxxxyxxy +−−+≥−++−+
49
Estas condições restringem as opções de valores que x
1
e x
2
podem assumir, visto
que, por exemplo, x
1
e x
2
não podem assumir valores (0,0).
Para contornar esta dificuldade, assume-se a alteração do Conjunto C para:
C’ = {x∈ R
n
: 1≤ x
1
≤ 50 e 0≤ x
i
≤ 50}.
E a solução do problema foi x = (1,0) como esperado. A rotina do algoritmo que roda
com dupla minimização simultânea está no anexo A.
5.3 – Terceiro Problema teste
Este problema foi adaptado do artigo de Marcotte e Zhu (1995).
Min x
2
– y
x
∈
C, y
∈
C
s.a.
Min y
2
+ 4 (5.13)
y
∈
C
Cujo conjunto C é dado por: C={x∈ R
2
: -1 ≤ x
i
≤ 1}.
Este é o primeiro exemplo que apresenta um problema de programação matemática
em dois níveis. A proposta é reformular o segundo nível usando uma função gap,
transformando o problema de dois níveis em outro de apenas um nível, com
soluções equivalentes.
Em primeiro lugar, deve-se transformar o problema de otimização do segundo nível,
usando sua condição de otimalidade de primeira ordem, em um problema de
desigualdade variacional parametrizada PDV (F(x, ⋅), C(x)):
Encontrar um y
∈
C(x) que satisfaz 〈 2y, (z - y) 〉 ≥ 0, ∀ z ∈ C(x)= [-1, 1] (5.14)
O segundo passo é formular a função gap associada a esta desigualdade variacional
definida por:
(5.15)
1
5212
2121
≥+≥+ xxexx
2
)(
||||..5,0)(,2max),(
zyzyyyxG
xCz
−−−=
∈
α
50
A função gap assim definida foi proposta pela primeira vez por Fukushima (1992),
chamada de função gap regularizada. Esta será a função gap utilizada para abordar
os demais exemplos, assumindo α = 1.
(5.16)
(5.17)
Como já é conhecido, a função gap é maior ou igual a zero, sob as hipóteses de
monotonia da F(x,⋅) e domínio C(x) compacto, e o mínimo de G(x, y*) é 0 onde y* é o
ponto ótimo do PDV (F(x, ⋅), C(x)); conseqüentemente, é a solução do problema do
segundo nível.
Sendo assim, podemos resolver o problema de dois níveis deste exemplo usando a
proposta da família de funções de penalização apresentada na seção 3.4.3,
transformando-o em um problema de apenas um nível:
Min x
2
– y + µ. G(x, y) (5.18)
x
∈
C,y
∈
C
onde µ é um inteiro positivo e G(x, y) é da forma (5.16).
Entretanto, para minimizar (5.18) usando a sub-rotina fmincon seria bom que se
tivesse a definição analítica da função gap. Isto se obtém executando o cálculo
indicado, isto é, precisa-se resolver:
(5.19)
A minimização tem como variável z, e a variável y é fixa para fazer a minimização.
Todavia, como já foi falado anteriormente, não se deve fixar apenas um valor de y e
pegar o resultado de z, assumindo-o como ponto ótimo. É preciso testar com vários
valores de y para se chegar à conclusão do comportamento de z em função do valor
de y.
Neste problema, testamos para:
• y = -1 que resultou no ponto z* = 1;
2
)(
2max),(
2
2
)(
yz
yyxG
xCz
+
−=
∈
+
−−−=
∈
2
)(
2min),(
2
2
)(
yz
yyxG
xCz
+
−−
∈
2
)(
2min
2
2
yz
y
Cz
51
• y = -0.3 que resultou no ponto z* = 0.3;
• y = 0.7 que resultou no ponto z* = -0.7;
• y = 1 que resultou no ponto z* = -1.
Logo, o mínimo (5.19) ocorrerá no ponto z* = -y, resultando a função gap
G(x, y) = -2y
2
.
Conhecendo este resultado, pode-se reescrever (5.18):
(5.20)
Agora, é possível minimizar (5.20) usando a sub-rotina fmincon, com x e y variáveis
e o parâmetro de penalização µ, que dá como resultado: o ponto ótimo (0, 1/4µ) que
converge para (0,0) quando µ tende ao infinito. Pode-se ver o código deste algoritmo
no Apêndice A1.
5.4 – Quarto Problema teste
Este problema de programação matemática em dois níveis linear foi retirado de
Vicente, L. (1992) (apud Dantas, S. (1998)).
Min x - 2y
x, y
s.a.
min y (5.21)
y
onde o conjunto viável C para o problema é a região definida pelas desigualdades
seguintes:
x + y ≤ 1
– x –y ≤ 1
x –y ≤ 1
– x +y ≤ 1
A região C delimita a variação da variável x ao intervalo [-1, 1]. Assim, fixada a
variável x, a variação da variável y depende dela, tendo como resultado para o
22
,
2min yyx
CyCx
µ
+−
∈∈
52
problema de segundo nível que o conjunto viável C(x) do segundo nível está definido
a seguir:
quando x ≥ 0 então x -1 ≤ y ≤ 1 – x e
quando x < 0 então -1 – x ≤ y ≤ 1 + x.
Pode-se notar que apesar de x ser variável apenas do primeiro nível, ela interfere no
conjunto de restrições do segundo nível. Assim, por exemplo, se x assume valores:
• x = 0.5, teremos z variando no intervalo [-0.5, 0.5];
• x = 0.9, teremos z variando no intervalo [-0.1, 0.1];
• x = 0.1, teremos z variando no intervalo [-0.9, 0.9];
• x = -0.2, teremos z variando no intervalo [-0.8, 0.8];
• x = -0.4, teremos z variando no intervalo [-0.6, 0.6].
Primeiro passo: transformar o problema de segundo nível em um PDV (F(x, ⋅), C(x)):
Encontrar um y ∈ C(x) que satisfaz 〈 1, (z - y) 〉 ≥ 0, ∀ z ∈ C(x). (5.22)
Segundo passo: formular a função gap associada a este problema de desigualdade
variacional que pode ser definida por:
,com α = 1 (5.23)
(5.24)
Sendo assim, pode-se formular a família de problemas de penalidade para
aproximar a solução do problema de dois níveis, penalizando a função gap,
transformando-o em um problema de apenas um nível, buscando forçar que G(x, y)
atinja o valor 0:
Min x – 2y + µ. G(x, y) (5.25)
x
∈
C,y
∈
C(x)
onde µ é um inteiro positivo e G(x, y) é calculado segundo (5.24).
2
)(
||||..5,0)(max),( zyzyyxG
xCz
−−−=
∈
α
( )
yzyzyxG
xCz
−+−−=
∈
2
)(
2
1
min),(
53
Entretanto, para minimizar (5.25) usando fmincon é preciso que se tenha os
resultados da minimização interna da função gap, isto é, precisa-se resolver:
(5.26)
Para encontrar a definição explícita de G(x, y), calcula-se a minimização, sendo a
variável da minimização z para valores fixos de y. Todavia, como já foi feito
anteriormente, deve-se testar vários valores de y para se chegar à conclusão do
comportamento de z em função do valor de y.
Observar que o conjunto de resposta do seguidor assume valores y não positivos
para cada x fixo, isto é, as respostas do seguidor são
Quando x<0 ⇒ y = -1 – x e quando x ≥ 0 ⇒ y = x – 1.
Agora, a verificação do comportamento de z, com y assumindo diferentes valores
não positivos:
Seja x = 0.2. Logo -0.8 ≤ z ≤ 0.8.
• y = -0.9 que resultou em z* = -0.8;
• y = -0.5 que resultou em z* = -0.8;
• y = -0.1 que resultou em z* = -0.8;
• y = 0 que resultou em z* = -0.8.
Logo, conclui-se que para qualquer y ∈ C’, a variável z assume valores z* = - (1- |x|)
e a solução de 5.26 é:
Agora, pode-se reescrever a família de problemas penalizados em (5.25) da forma:
(5.28)
( )
yzyzyxG
xCz
−+−=
∈
2
)(
2
1
min),(
( ) ( )
−−+−+−=
−+−
∈
yxyxyzyz
Cz
1||||1
2
1
2
1
min
22
( )
−−+−+−−+−
∈∈
yxyxyx
CyCx
1||||1
2
1
.2min
2
,
µ
54
E minimizar (5.28) usando fmincon, com x e y variáveis, encontrando o ponto ótimo
(x*,y*) = (-1, 0). O código deste algoritmo encontra-se no Apêndice A2. A não
diferenciabilidade na formulação (5.28) se evita rodando dois problemas, finalmente
escolhendo o menor.
5.4 – O Problema das redes de transporte de gás natural
O modelo do problema das redes de transporte de gás natural, proposto por
Kalashnikov e Ríos-Mercado (2006), foi detalhado no capítulo 2 e pré-processado no
capítulo 4.
Este problema é primeiramente modelado como um problema de programação
matemática em dois níveis, que pode ser escrito da forma:
min h
1
(x,s,y,u,v,z) = z
s.a.
(2.1b)-(2.1e) (5.29)
min h
3
(x,s,y,u,v,z) = |z|
s.a.
(2.2b)-(2.2l)
Entretanto, para apropriá-lo aos métodos de solução existentes para modelos de
programação matemática de apenas um nível, este problema foi reformulado no
capítulo 4, onde foi feito o pré-processamento, mantendo as equivalências em
relação às soluções.
Para que se tenha uma melhor compreensão da maneira como o algoritmo foi
montado, apresenta-se a seguir o diagrama de blocos deste modelo:
55
Regularizaçao
Regularização da função do 2
°
nível (|z|
→
z
2
)
Alteração da função objetivo (z
→
F = -z)
Retirada da restrição
2.2 (l). Qual é o valor
de q?
O valor de q é: q=0
Transformação do problema
do 2
° nível em um PDV (F,C)
C
C)
Transformação do PDV (F,C)
em um problema de
otimização equivalente
usando uma função gap
Uso da função gap como
termo de penalidade para o
problema do 1
° nível
O valor de q é: q=1
Transformação do problema
do 2
° nível em um PDV (F,C)
C)
Transformação do PDV (F,C)
em um problema de
otimização equivalente
usando uma função gap
Uso da função gap como
termo de penalidade para o
problema do 1
° nível
Comparação dos
resultados. Qual é o
menor valor?
O menor valor é a solução do problema
original
56
O problema pré- processado, apresentado em (4.4), pode ser reescrito da forma:
(5.30)
onde µ é um número positivo, α é um número não-negativo e o subconjunto C
β
é
dado por: C
β
= {( w
1
,w
2
β
) : (2.1b) – (2.1e), (2.2b) – (2.2h), (2.2j) – (2.2k), q = β}.
Para minimizar (5.30) Kalashnikov e Ríos-Mercado sugerem um exemplo de um
pequeno problema teste deste modelo para 4 zonas e dois dias, isto é, para N = 2 e
P = 4. Com os seguintes dados de entrada:
Tabela 1: Desequilíbrios iniciais e coeficientes
Zonas Desequilíbrio inicial x
0
Fator de penalidade r
j
Coeficiente δ
j
1 -10 10 0,1
2 -4 8 0,2
3 3 6 0,3
4 6 4 0,4
Tabela 2: Limites superiores e inferiores para o total de desequilíbrios diários
Dias Limite inferior (x
t
l
) Limite superior ( x
t
u
)
1 -11 1
2 -17 7
Tabela 3: Limites inferiores para cada x
ti
Dias Zona 1 Zona 2 Zona 3 Zona 4
1 -15 -9 0 1
2 -20 -13 -5 -4
Tabela 4: Limites superiores para cada x
ti
Dias Zona 1 Zona 2 Zona 3 Zona 4
1 -5 1 6 10
2 0 5 10 15
−+−∇−+−
∈∈
,
2
1
),()(min.)(min
2
222222
2
2
,
22
1
βββββ
αµ
β
β
β
wwwwwFwFwF
WwCww
57
Tabela 5: Custos por transportar (empurrar) gás entre zonas (f
ij
)
Zona 1 Zona 2 Zona 3 Zona 4
Zona 1 2 4 4
Zona 2 2 2
Zona 3 1
Zona 4
Tabela 6: Porcentagem de gás retido como combustível quando transportado entre
zonas (e
ij
)
Zona 1 Zona 2 Zona 3 Zona 4
Zona 1 0,1 0,2 0,2
Zona 2 0,1 0,1
Zona 3 0,05
Zona 4
Tabela 7: Bonificação do emissor por puxar gás entre zonas (b
ij
)
Zona 1 Zona 2 Zona 3 Zona 4
Zona 1 4 2 2
Zona 2 4 4
Zona 3 2
Zona 4
Tabela 8: Limites inferiores para cada s
ti
Dias Zona 1 Zona 2 Zona 3 Zona 4
1 -3 -3 -3 -3
2 -3 -3 -3 -3
Tabela 9: Limites superiores para cada s
ti
Dias Zona 1 Zona 2 Zona 3 Zona 4
1 3 3 3 3
2 3 3 3 3
58
Kalashnikov e Ríos-Mercado montaram um algoritmo de busca direta para verificar a
eficiência do pré-processamento, no qual o problema de dois níveis, cujo modelo foi
detalhado no capítulo 2, transformou-se em um problema de um nível com
penalização inexata.
Os resultados obtidos para x
ti
e s
ti
ótimos apresentam-se nas tabelas 10 e 11:
Tabela 10: Valores ótimos para as variações de desequilíbrio (s
ti
)
Dias Zona 1 Zona 2 Zona 3 Zona 4
1 3 -3 3 3
2 2.999 -2,998 3 2,999
Tabela 11: Valores ótimos para os desequilíbrios diários (x
ti
)
Dias Zona 1 Zona 2 Zona 3 Zona 4
1 -7 -7 6 9
2 -4,001 -9,998 9 11,999
Com y = (0,0,0,7), u = (0,0,0,0,0,0), v = (0,0,4.001,9,0.998,0) e o valor ótimo da
função objetivo z = -75.994.
Usando fmincon, foi possível chegar ao seguinte resultado:
Tabela 12: Valores ótimos para as variações de desequilíbrio (s
ti
)
Dias Zona 1 Zona 2 Zona 3 Zona 4
1 3 -2,3741 1,3741 3
2 3 -3 -3 3
Tabela 13: Valores ótimos para os desequilíbrios diários (x
ti
)
Dias Zona 1 Zona 2 Zona 3 Zona 4
1 -7 -6,3740 4,3740 9
2 -4,001 -9,3740 1,3740 12
y = (0, 0, 0, 0), u = (0,0,0,0,0,0), v = (0, 1.3741, 2.6259, 0, 9.3741, 0) e o valor ótimo
da função objetivo z = - 45.4963.
59
Ao comparar estas respostas com os resultados numéricos do trabalho de
Kalashnikov e Ríos-Mercado (2006), é possível perceber que os valores das
variáveis x
ti
e s
ti
são próximos, exceto na zona 3. Os desequilíbrios diários na zona 3
alteram o desequilíbrio final da zona 4 e, conseqüentemente, os fluxos transladados
entre zonas.
60
Capítulo 6
CONCLUSÕES
Neste trabalho foi abordado o problema de minimizar custos penalizados na
distribuição de gás natural. Inicialmente ele é modelado como um problema de
programação matemática não linear em dois níveis misto-inteiro. Pelo fato de se ter
apenas uma variável booleana, opta-se por resolver dois problemas assumindo os
valores 0 e 1 para esta variável. Para facilitar a sua resolução, este modelo foi pré-
processado, usando desigualdades variacionais, funções gap e penalidades, o que o
transforma em problema de programação matemática de apenas um nível, usando a
função gap como um termo de penalidade para o problema de primeiro nível.
Os conceitos de desigualdades variacionais, funções gap, métodos de penalidades e
problemas de programação de dois níveis foram estudados a fim de entender os
resultados teóricos necessários para a abordagem deste modelo do problema de
distribuição de gás natural focado. Este estudo foi muito útil para o pré-
processamento do modelo e para a implementação do algoritmo proposto.
Para verificar o método aplicado, foram usados problemas testes que, apesar de
dimensões bem menores que as dimensões do problema das redes de gás natural,
foram de grande valia para o aprendizado e aplicação das técnicas de pré-
processamento sugeridas no capítulo 4. Com eles foi possível implementar
algoritmos, usando o ambiente MATLAB, e verificar o comportamento da função gap.
O exemplo do problema das redes de transporte de gás natural, para 4 zonas e 2
dias, sugerido por Kalashnikov e Ríos-Mercado também foi implementado usando o
ambiente MATLAB. O objetivo de rodar este exemplo com o método proposto foi a
verificação dos resultados relatados no artigo, porém, ainda fica para trabalhos
futuros outros testes a serem rodados que fornecerão maiores informações para
comparar com os resultados numéricos do artigo trabalhado.
61
O desenvolvimento deste trabalho serviu para aplicar e contextualizar conceitos mais
elaborados da matemática e de cálculo numérico para resolver problemas concretos
do planejamento da produção e exploração de petróleo e derivados, assim como
também se conseguiu o aprofundamento no uso do ambiente MATLAB nas rotinas
de otimização.
Como trabalhos futuros, sugere-se a pesquisa de algoritmos de penalidades exatas
para desigualdades variacionais associadas à função gap inicialmente sugerida por
Hearn (1981). Adicionalmente sugere-se tentar o uso de algoritmos meta-heurísticos
na solução de problema de dois níveis, usando discretizações nas variáveis, o que
aumentaria o tamanho do problema, porém, em compensação, essas heurísticas
liberam da condição de diferenciabilidade sobre as funções.
62
7 – Referências Bibliográficas
1. André, T. A. (2007) Penalidades Exatas para Desigualdades Variacionais.
Tese (Mestrado em Ciências) – São Paulo, SP – Instituto de Matemática e
Estatística da Universidade de São Paulo – USP, 67p.
2. Aubin, J. P. (1998) Optima and Equilibria An Introduction to Nonlinear
Analysis. 2., Springer.
3. Auchmuty, G. (1989) Variational Principles for Variational Inequalities.
Numerical Functional Analysis and Optimization, Texas, 10:863-874.
4. Auslender, A.; (1976) Optimization: méthodes numériques. Masson S. A. 120
Bd Saint-Germain, Paris 6.
5. Bazaraa, M. S.; Sherali, H. D.; Shetty, C. M. (1993) Nonlinear Programming
Theory and Algorithms. USA: A Wiley- Interscience publication, New York,
EUA.
6. Berkovitz, L. D. (2001) Convexity and Optimization in R
n
. USA : A Wiley-
Interscience publication.
7. Dafermos, S. C. (1980), Treffic Equilibrium and Variational inequalities.
Transportation Science, 14:42-54.
8. Dantas, S. (1998) Problemas de Programação de Dois Níveis: um estudo dos
casos Linear, Linear-quadrático e Quadrático. Dissertação de Mestrado em
Ciências em engenharia de Sistemas de Computação. COPPE/UFRJ, Rio de
Janeiro.
9. Dempe, S.; Kalashnikov, V.; Ríos-Mercado R. Z. (2005) Discrete bilevel
programing: Application to a natural gas cash-out problem. European Journal
of Operation Research 166, 469-488.
10. Fukushima M. (1992) Equivalent differentiable optimization problems and
descents methods for asymmetric variational inequality problems.
Mathematical Programming 53, 99-110.
11. Gale, D.; (1967) A Geometric Duality Theorem with Economic Applications.
Review of Economic Studies, 34, pp. 19-24.
12. Hartman, P; Stampacchia, G. (1966) On Some Nonlinear Elliptic Differential
Functional Equations. Acta Mathematica, 115:153-188.
63
13. Hearn, D. W. (1981) The Gap Function of a Convex Program. Operations
Research Letters, USA, 67-71.
14. Hiriart-Urruty, J B.; Lemaréchal, C. (1993) Convex Analysis and Minimization
Algorithms Vol.I . Springer-Verlag.
15. Iamashita, E. K.; Galaxe, F.; Arica, J.; Justiniano, L. R. S.; Iachan, R. (2005)
Um Algoritmo Genético Híbrido para o Planejamento de Movimentação de
Gás da Bacia de Campos. In: Anais do XXXVII SBPO - Simpósio Brasileiro de
Pesquisa Operacional, GRAMADO. v. 1. p. 1-12.
16. Izmailov, A.; Solodov, M. (2005) Otimização - volume 1: Condições de
Otimalidade, Elementos de Análise Convexa e de Dualidade. Rio de Janeiro:
Impa, 253p.
17. Kalashnikov, V. V.; Ríos-Mercado, R. Z. (2006) A Natural Gás Cash-Out
Problem: A Bilevel Programming Framework and a Penalty Function Method.
México, 1-19.
18. Larsson, T.; Patriksson, M. (1994) A class of gap functions for variational
inequalities. Mathematical Programming 64:53-79.
19. Marcotte, P.; Zhu, D. L. (1996) Exact and Inexact Penalty Methods for the
Generalized Bilevel Programming Problem. Mathematical Programming, 74:
141-157.
20. Marcotte P.; Dussault J. P. (1989) A sequential linear programming algorithm
for solving monotone variational inequalities. SIAM Journal of Control and
Optimization 27, 1260-1278.
21. Morales, G. (1997) O problema de programação matemática com restrições
generalizadas de Equilíbrio. Tese (Doutorado em Engenharia de Sistemas de
Computação) – Rio de Janeiro, RJ – Universidade Federal do Rio de Janeiro -
Coppe/UFRJ.
22. Morales, Y. V.: Ríos-Mercado, R. Z. (2002) Preprocessamiento Efectivo de un
Problema de Minimización de Combustible en Sistemas de Transporte de Gas
Natural. División de Posgrado em Ingeniería de Sistemas – Serie de Reportes
Técnicos, México, 1-19.
23.
Ozdaglar A. E., Pseudonormality and a Lagrange Multiplier Theory for
Constrained Optimization. Ph.D. Thesis, January 2003, MIT, EECS
Department.
64
24. Ríos-Mercado, R. Z.; Kim, S.; Boyd, E. A. (2006) Efficient Operation of Natural
Gas Transmission Systems: A Network- Based Heurístic for Cyclic Structures.
Computers & Operations Research, 23(8):2323-2351.
25. Ríos-Mercado, R. Z.; Wu, S.; Scott, L. R.; Boyd, E. A. (2000) Preprocessing
on Natural Gas. Transmission Networks, Reporte técnico, UANL, San Nicolás
de los Garza, México.
26. Smith, M. J. (1979) The existence, uniqueness and stability of traffic equilibria.
Transportation Reserach, 13:295-304.
27. Souza, J. S. (2005) Alguns tópicos teóricos e algorítmicos relacionados a um
modelo de Programação Linear em Dois Níveis. Tese (Mestrado em
Engenharia de Produção) – Campos dos Goytacazes – RJ, Universidade
Estadual do Norte Fluminense - UENF, 13-31p.
28. Trujillo, C. E. V. (2003) Algoritmos de Descida baseados em Funções Gap
para resolver Problemas Variacionais Monótonos. Tese (Mestrado em
Engenharia de Produção) – Campos dos Goytacazes – RJ, Universidade
Estadual do Norte Fluminense - UENF, 67p.
29.
Universidade Estadual do Norte Fluminense - Câmara de Pesquisa e Pós-
Graduação (1997) Normas para a Elaboração de Tese na Uenf: Resolução
n° 006/97. Campos dos Goytacazes.
30. Vicente, L. N. (1992) Programação de Dois Níveis. Tese de Mestrado,
Departamento de Matemática, Universidade de Coimbra, Portugal.
31. Wu, S.; Ríos-Mercado, R. Z.; Boyd, E. A.; Scott, L. R. (2000) Model
Relaxations for the Fuel Cost Minimization of Steady-State Gas Pipeline
Networks. Mathematical and Computer Modelling, 31(2-3): 197-220.
32. Zuhovickii, S. I.; Poljak R. A.; Primak M. E. (1969a) Two methods of search for
equilibrium points of n-person concave games. Soviet Mathematics Doklady
10, 279-282.
33. Zuhovickii, S. I.; Poljak R. A.; Primak M. E. (1969b) Numerical Methods of
finding equilibrium points of n-person games. Proceedings of the First Winter
School of mathematical Programming at Drogobych 1, 93-130.
34. Zuhovickii, S. I.; Poljak R. A.; Primak M. E. (1970) On a n-person concave
game and a production model. Soviet Mathematics Doklady 11, 522-526.
35. Zuhovickii, S. I.; Poljak R. A.; Primak M. E. (1973) Concave multiperson
games: numerical methods. Matekon 9, 10-30.
65
Anexo A
Para o desenvolvimento deste trabalho foram necessários alguns conceitos básicos
e resultados de análise convexa colocados a seguir.
A.1 - Função Convexa
Seja C um conjunto não vazio e convexo de R
n
e f: C→R uma função dada. A
função f é convexa sobre C quando:
f(αx + (1-α) x’) ≤ α f(x) + (1-α) f(x’)
sempre que x, x’ ∈ C, α ∈ (0,1).
A.2 - Função Côncava
Uma função f: C→R é côncava se (– f) é uma função convexa em C.
A.3 - Função Inferiormente Semicontínua
Seja C um subconjunto de R
n
. A função f:C→R é inferiormente semicontínua no
ponto x ∈ C quando para qualquer seqüência {x
k
} ⊂ C tal que {x
k
}→x (k→∞), tem-se:
A função f é inferiormente semicontínua em C, se é inferiormente semicontínua em
cada ponto x ∈ C.
A.4 - Função Conjugada
Em particular, dada a função estendida f: R
n
→ R ∪ {∞}, sua função conjugada é
f*: R
n
→ R ∪ {∞}, com:
A.5 - Domínio Essencial
Quando f: R
n
→ R ∪ {∞} é uma função, seu domínio essencial é:
{
}
)(,sup)(* xfyxyf
n
Rx
−=
∈
)()(inflim xfxf
k
k
≥
∞→
66
dom f = {x ∈ R
n
: |f(x)| < ∞ }.
A.6 - Função Indicadora
Seja C um subconjunto de R
n
. A função indicadora de C é a função I
C
: R
n
→R ∪ {∞}
definida por:
A.7 - Função Suporte
A função conjugada da função indicadora I
C
é chamada de função suporte de C,
denotada por:
A.8 – Subdiferencial
Uma relaxação da condição de derivabilidade de uma função convexa, definida
sobre R
n
, leva a definir subgradiente.
Seja f: R
n
→R uma função convexa. Dizemos que y ∈ R
n
é um subgradiente de f no
ponto x ∈ R
n
se
f(z) ≥ f(x) + 〈z – x,y〉 ; ∀ z ∈ R
n
.
O conjunto de todos os subgradientes de f em x diz-se o subdiferencial de f em x;
denotado por ∂f(x).
A.9 - Função Localmente Lipschitz
Uma função f: R
n
→R ∪ {∞} é localmente Lipschitz em um subconjunto aberto Ω ⊂ R
n
se para cada x
0
∈ Ω, existe um η > 0 e um c > 0 tal que
∀ x,y ∈ B(x
0
,η), |f(x) – f(y)| ≤ c ||x – y||
∞
∈
=
contráriocaso
Cxse
xI
C
0
)(
[
]
)(,sup)(
*
xIyxyI
C
Cx
C
−=
∈
yxyI
Cx
C
,sup)(
*
∈
=
67
A.10 - Generalização da Desigualdade de Young
Os princípios variacionais utilizados baseiam-se na generalização da desigualdade
de Young definida em [12].
Suponha f: R
n
→ R ∪ {∞} uma função convexa inferiormente semicontínua e f* sua
função conjugada, então:
f(x) + f*(y) ≥ 〈x,y〉 para todo x,y em R
n
A igualdade acontece se, e somente se y ∈∂ f (x).
A.11 - Generalização da Desigualdade de Cauchy
Quando C é um subconjunto de R
n
não vazio e limitado, usando a desigualdade de
Cauchy |〈z,w〉| ≤ |z| |w| e a relação (4), temos:
|〈x,y〉| ≤ ||x|| ||y||
Se ||x|| < k ∀ x ∈ C e k > 0
|〈x,y〉| ≤ k ||y||
||||,sup ykyx
Cx
≤
∈
|I*
C
(y)| ≤ k ||y||
A.12 - Cone Convexo
• Seja C um subconjunto não-vazio e convexo de R
n
. C é chamado cone
convexo se: qualquer que seja x ∈ C ⇒λx ∈ C para todo λ≥0.
• O cone conjugado do conjunto C se define por:
C* = {y ∈ R
n
: 〈x,y〉 ≥ 0 para todo x ∈ C}.
{
}
yxxfyxxf
n
Rx
,)(,sup)( ≥−+
∈
{
}
0)](,[)(,sup ≥−−−
∈
xfyxxfyx
n
Rx
68
Então, temos que o conjunto:
(i) (- C*) = {y ∈ R
n
: 〈x,y〉 ≤ 0 para todo x ∈ C} conhecido como cone polar ou
também como cone dual [14].
(ii)
(iii)
Da relação (4), pelo fato do conjunto C ser um cone convexo e fechado, pode-se
concluir que:
I*
C
(y) = I
-C*
(y)
No caso em que é necessário construir o cone relativo a um ponto na fronteira (ou
fecho) de um conjunto convexo C, define-se o cone normal.
• O cone normal do conjunto no ponto x* se define como:
N
C
(x*)= {v∈ R
n
: 〈v, y-x*〉 ≤ 0 para todo y ∈ C}
∞
−∈
=
−
contráriocaso
Cyse
yI
C
*0
)(
*
yxyI
Cx
C
,sup)(
*
∈
−
=
69
Apêndice A
A1 - Programa principal: exemplo_marcotte.m
global lb;
global ub;
global m;
lb = [-1;-1];
ub = [1;1];
m = 100;
x0 = [1;0];
[x,resul,EXITFLAG,OUTPUT] = fmincon(@funcaomarcote,x0,[],[],[],[],lb,ub)
Sub-rotina_A1: funcaomarcote.m
function fm=funcaomarcote(x)
global m;
fm = x(1)^2 - x(2)+ m*2*x(2)^2;
A2 – Programa principal: exemplo_dantas.m
global A;
global B;
global m;
A = [1 1;-1 -1;1 -1;-1 1];
B = [1;1;1;1];
m = 10;
w0 = [-0.2;-0.8];
[w,resul,EXITFLAG,OUTPUT] = fmincon(@funcaodantaspenalizado,w0,A,B)
Sub-rotina_A2: funcaodantaspenalizado.m
function fdp=funcaodantaspenalizado(w)
global A;
global B;
global m;
fs = funcaodantas(w)
fdg = funcaodantasgap(w)
fsp = fs + m* fsg
Sub-rotina_A2: funcaodantas.m
function fd = funcaodantas(w)
fd = w(1)-2*w(2);
Sub-rotina_A2: funcaodantasgap.m
function fdg=funcaodantasgap(w)
z = -(1-abs(w(1)))
fdg = -(0.5*(z-w(2))^2 + z-w(2));
70
A3 – Exemplo principal: exemplo_roger_slava_penalizado.m
global a
global M
global q
global r
global d
global e
global b
global f
global A
global B
global Aeq
global Beq
global snew
snew = [-8 -1 0 5 -11 2 -2 5];
xt3 = [-10 -1 0 4 -8 2 1 2];
st3 = [0 3 -3 -2 2 3 1 -2];
y3 = [-3 0 0 0];
u3 = [0 0 0 0 0 0];
v3 = [2 1 2 0 0 0];
xt4 = [-8 -1 0 5 -11 2 -2 5];
st4 = [2 3 -3 -1 -3 3 -2 0];
y4 = [-6 0 -0.1 0];
u4 = [0 0 0 1 0 0];
v4 = [1 0 4 0 0 1];
w2b =[xt3 st3 y3 u3 v3];
w0 = [xt4 st4 y4 u4 v4];
a = 9;
M = 100;
q = 0;
e = [0.1 0.2 0.2 0.1 0.1 0.05];
%d = [0.1 0.2 0.3 0.4];
d = [0 0 0 0];
r = [10 8 6 4];
b = [4 2 2 4 4 2];
f = [2 4 4 2 2 1];
Mq = M*q;
% Agora define-se A1n (56x16) como a matriz 0
A1n = 0
%restricao 1b
A1n(1,:) = [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
A1n(2,:) = [0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
A1n(3,:) = [0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0];
A1n(4,:) = [0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0];
A1n(5,:) = [0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0];
A1n(6,:) = [0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0];
71
A1n(7,:) = [0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0];
A1n(8,:) = [0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0];
A1n(9,:) = [-1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
A1n(10,:) = [0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
A1n(11,:) = [0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0];
A1n(12,:) = [0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0];
A1n(13,:) = [0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0];
A1n(14,:) = [0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0];
A1n(15,:) = [0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0];
A1n(16,:) = [0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0];
%restricao 1c
A1n(17,:) = [0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0];
A1n(18,:) = [0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0];
A1n(19,:) = [0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0];
A1n(20,:) = [0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0];
A1n(21,:) = [0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0];
A1n(22,:) = [0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0];
A1n(23,:) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0];
A1n(24,:) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1];
A1n(25,:) = [0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0];
A1n(26,:) = [0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0];
A1n(27,:) = [0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0];
A1n(28,:) = [0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0];
A1n(29,:) = [0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0];
A1n(30,:) = [0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0];
A1n(31,:) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0];
A1n(32,:) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1];
%restricao 1d
A1n(33,:) = [1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0];
A1n(34,:) = [0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0];
A1n(35,:) = [-1 -1 -1 -1 0 0 0 0 0 0 0 0 0 0 0 0];
A1n(36,:) = [0 0 0 0 -1 -1 -1 -1 0 0 0 0 0 0 0 0];
%Agora montar a matriz A2n 56x16 como a matriz 0
A2n = 0
%restricao 2h
A2n(37,:) = [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
A2n(38,:) = [0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
A2n(39,:) = [0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0];
A2n(40,:) = [0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0];
A2n(41,:) = [-1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
A2n(42,:) = [0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
A2n(43,:) = [0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0];
A2n(44,:) = [0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0];
%restricao 2k
A2n(45,:) = [0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0];
A2n(46,:) = [0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0];
A2n(47,:) = [0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0];
A2n(48,:) = [0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0];
A2n(49,:) = [0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0];
A2n(50,:) = [0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0];
72
A2n(51,:) = [0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0];
A2n(52,:) = [0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0];
A2n(53,:) = [0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0];
A2n(54,:) = [0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0];
A2n(55,:) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0];
A2n(56,:) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1];
A = [A1n A2n];
% Agora definiremos o vetor B(56x1)
%B =
[0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;Mq;Mq;Mq;Mq;M-
Mq;M-Mq;M-Mq;M-Mq;0;0;0;0;0;0;0;0;0;0;0;0];
B = [-5;1;6;10;0;5;10;15;15;9;0;-
1;20;13;5;4;3;3;3;3;3;3;3;3;3;3;3;3;3;3;3;3;1;7;11;17;Mq;Mq;Mq;Mq;M-Mq;M-Mq;M-
Mq;M-Mq;0;0;0;0;0;0;0;0;0;0;0;0];
% Agora vamos definir a matriz Aeq1n(13x16) como a matriz 0
Aeq1n=0
Aeq1n(1,:)= [0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0];
Aeq1n(2,:)= [0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0];
Aeq1n(3,:)= [0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0];
Aeq1n(4,:)= [0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0];
Aeq1n(5,:)= [0 0 0 0 -1 -1 -1 -1 0 0 0 0 0 0 0 0];
Aeq1n(6,:)= [1 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0];
Aeq1n(7,:)= [0 1 0 0 0 0 0 0 0 -1 0 0 0 0 0 0];
Aeq1n(8,:)= [0 0 1 0 0 0 0 0 0 0 -1 0 0 0 0 0];
Aeq1n(9,:)= [0 0 0 1 0 0 0 0 0 0 0 -1 0 0 0 0];
Aeq1n(10,:)= [0 0 0 0 1 0 0 0 -1 0 0 0 -1 0 0 0];
Aeq1n(11,:)= [0 0 0 0 0 1 0 0 0 -1 0 0 0 -1 0 0];
Aeq1n(12,:)= [0 0 0 0 0 0 1 0 0 0 -1 0 0 0 -1 0];
Aeq1n(13,:)= [0 0 0 0 0 0 0 1 0 0 0 -1 0 0 0 -1];
% Agora definir a matriz Aeq2n(13x16) como a matriz 0
Aeq2n=0
Aeq2n(1,:)= [1 0 0 0 1 1 1 0 0 0 -1 -1 -1 0 0 0];
Aeq2n(2,:)= [0 1 0 0 (e(1,1)-1) 0 0 1 1 0 1 0 0 -1 -1 0];
Aeq2n(3,:)= [0 0 1 0 0 (e(1,2)-1) 0 (e(1,4)-1) 0 1 0 1 0 1 0 -1];
Aeq2n(4,:)= [0 0 0 1 0 0 (e(1,3)-1) 0 (e(1,5)-1) (e(1,6)-1) 0 0 1 0 1 1];
Aeq2n(5,:)= [1 1 1 1 e(1,1) e(1,2) e(1,3) e(1,4) e(1,5) e(1,6) 0 0 0 0 0 0];
Aeq = [ Aeq1n Aeq2n];
% Agora definir o vetor Beq(13x1)
Beq = [0;0;0;0;0;-10;-4;3;6;-10;-4;3;6];
[w2b,fvalpenal,EXITFLAG,OUTPUT]=fmincon(@funcaopenalizada_roger_slava,w0,A,B,Aeq,Beq,[],[],@nonlcon)
Sub-rotina_A3: funcaopenalizada_roger_slava.m
global a
73
global M
global q
global r
global d
global e
global b
global f
global A
global B
global Aeq
global Beq
global snew
xt3 = [-10 -1 0 4 -8 2 1 2];
st3 = [0 3 -3 -2 2 3 1 -2];
y3 = [-3 0 0 0];
u3 = [0 0 0 0 0 0];
v3 = [2 1 2 0 0 0];
xt4 = [-8 -1 0 5 -11 2 -2 5];
st4 = [2 3 -3 -1 -3 3 -2 0];
y4 = [-6 0 -0.1 0];
u4 = [0 0 0 1 0 0];
v4 = [1 0 4 0 0 1];
Mq = M*q;
w2bnew = w2b(1,17:32);
ff = funcaof_roger_slava (w2b)
fgnew = funcaogapnew (w2bnew)
w2b(1,17:32)= w2bnew;
fp = ff + M*abs(fgnew)
snew = w2b(1,1:8);
Sub-rotina_A3: funcaof_roger_slava.m
function ff = funcaof_roger_slava(w2b)
global a
global M
global q
global r
global d
global e
global b
global f
global A
global B
global Aeq
global Beq
ff = -(r*[w2b(1,17) w2b(1,18) w2b(1,19) w2b(1,20)]' - d*[(max([0 w2b(1,17)]))^2 (max([0
w2b(1,18)]))^2 (max([0 w2b(1,19)]))^2 (max([0 w2b(1,20)]))^2]' + b*[w2b(1,27) w2b(1,28)
w2b(1,29) w2b(1,30) w2b(1,31) w2b(1,32)]' - f*[w2b(1,21) w2b(1,22) w2b(1,23) w2b(1,24)
74
w2b(1,25) w2b(1,26)]'+ f(1,1)*e(1,1)*w2b(1,21) + f(1,2)*e(1,2)*w2b(1,22) +
f(1,3)*e(1,3)*w2b(1,23) + f(1,4)*e(1,4)*w2b(1,24) + f(1,5)*e(1,5)*w2b(1,25) +
f(1,6)*e(1,6)*w2b(1,26));
Sub-rotina_A3: funcaogapnew.m
function fgnew = funcaogapnew(w2bnew)
global a
global M
global q
global r
global d
global e
global b
global f
global snew
xt3 = [-10 -1 0 4 -8 2 1 2];
st3 = [0 3 -3 -2 2 3 1 -2];
y3 = [-3 0 0 0];
u3 = [0 0 0 0 0 0];
v3 = [2 1 2 0 0 0];
xt4 = [-8 -1 0 5 -11 2 -2 5];
st4 = [2 3 -3 -1 -3 3 -2 0];
y4 = [-6 0 -0.1 0];
u4 = [0 0 0 1 0 0];
v4 = [1 0 4 0 0 1];
w02n = [y4 u4 v4];
Mq = M*q;
%Agora se monta a matriz A2n 20x16 (restricoes 2h e 2k)
A2n =[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 ;
-1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 ;
0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 ;%inicio da restricao 2k
0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 ;
0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 ;
0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 ;
0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 ;
0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 ;
0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 ;
0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 ;
0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 ;
0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 ;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 ;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 ];
75
% Agora se define o vetor B(20x1)
B2n = [Mq;Mq;Mq;Mq;M-Mq;M-Mq;M-Mq;M-Mq;0;0;0;0;0;0;0;0;0;0;0;0];
% Agora definir a matriz Aeq2n(5x16)
Aeq2n = [1 0 0 0 1 1 1 0 0 0 -1 -1 -1 0 0 0;
0 1 0 0 (e(1,1)-1) 0 0 1 1 0 1 0 0 -1 -1 0;
0 0 1 0 0 (e(1,2)-1) 0 (e(1,4)-1) 0 1 0 1 0 1 0 -1;
0 0 0 1 0 0 (e(1,3)-1) 0 (e(1,5)-1) (e(1,6)-1) 0 0 1 0 1 1;
1 1 1 1 e(1,1) e(1,2) e(1,3) e(1,4) e(1,5) e(1,6) 0 0 0 0 0 0];
% Agora vamos definir o vetor Beq(5x1)
Beq2n=[snew(1,5);snew(1,6);snew(1,7);snew(1,8);snew(1,5)+snew(1,6)+snew(1,7)+snew(1,8
)];
[w2new,fvalinicial,EXITFLAG,OUTPUT]=fmincon(@funcaogap_roger_slavanew,w02n,A2n
,B2n,Aeq2n,Beq2n,[],[],@nonlconnew)
fgnew = -(((r*[w2bnew(1,1) w2bnew(1,2) w2bnew(1,3) w2bnew(1,4)]' - d*[(max([0
w2bnew(1,1)]))^2 (max([0 w2bnew(1,2)]))^2 (max([0 w2bnew(1,3)]))^2 (max([0
w2bnew(1,4)]))^2]' + b*[w2bnew(1,11) w2bnew(1,12) w2bnew(1,13) w2bnew(1,14)
w2bnew(1,15) w2bnew(1,16)]' - f*[w2bnew(1,5) w2bnew(1,6) w2bnew(1,7) w2bnew(1,8)
w2bnew(1,9) w2bnew(1,10)]'+ f(1,1)*e(1,1)*w2bnew(1,5) + f(1,2)*e(1,2)*w2bnew(1,6) +
f(1,3)*e(1,3)*w2bnew(1,7) + f(1,4)*e(1,4)*w2bnew(1,8) + f(1,5)*e(1,5)*w2bnew(1,9) +
f(1,6)*e(1,6)*w2bnew(1,10))*[ r(1,1)-2*d(1,1)*max([0 w2bnew(1,1)]) r(1,2)-
2*d(1,2)*max([0 w2bnew(1,2)]) r(1,3)-2*d(1,3)*max([0 w2bnew(1,3)]) r(1,4)-
2*d(1,4)*max([0 w2bnew(1,4)]) -f(1,1)*(1 - e(1,1)) -f(1,2)*(1 - e(1,2)) -f(1,3)*(1 - e(1,3)) -
f(1,4)*(1 - e(1,4)) -f(1,5)*(1 - e(1,5)) -f(1,6)*(1 - e(1,6)) b(1,1) b(1,2) b(1,3) b(1,4) b(1,5)
b(1,6)])*(w2new - w2bnew)' + 0.5* a * (norm(w2new - w2bnew))^2);
Sub-rotina_A3: funcaogap_roger_slavanew
function fgn = funcaogap_roger_slavanew(w2new)
global a
global M
global q
global r
global d
global e
global b
global f
global snew
xt3 = [-10 -1 0 4 -8 2 1 2];
st3 = [0 3 -3 -2 2 3 1 -2];
y3 = [-3 0 0 0];
u3 = [0 0 0 0 0 0];
v3 = [2 1 2 0 0 0];
xt4 = [-8 -1 0 5 -11 2 -2 5];
st4 = [2 3 -3 -1 -3 3 -2 0];
y4 = [-6 0 -0.1 0];
u4 = [0 0 0 1 0 0];
76
v4 = [1 0 4 0 0 1];
w2bnew = [y3 u3 v3];
fgn = ((r*[w2bnew(1,1) w2bnew(1,2) w2bnew(1,3) w2bnew(1,4)]' - d*[(max([0
w2bnew(1,1)]))^2 (max([0 w2bnew(1,2)]))^2 (max([0 w2bnew(1,3)]))^2 (max([0
w2bnew(1,4)]))^2]' + b*[w2bnew(1,11) w2bnew(1,12) w2bnew(1,13) w2bnew(1,14)
w2bnew(1,15) w2bnew(1,16)]' - f*[w2bnew(1,5) w2bnew(1,6) w2bnew(1,7) w2bnew(1,8)
w2bnew(1,9) w2bnew(1,10)]'+ f(1,1)*e(1,1)*w2bnew(1,5) + f(1,2)*e(1,2)*w2bnew(1,6) +
f(1,3)*e(1,3)*w2bnew(1,7) + f(1,4)*e(1,4)*w2bnew(1,8) + f(1,5)*e(1,5)*w2bnew(1,9) +
f(1,6)*e(1,6)*w2bnew(1,10))*[ r(1,1)-2*d(1,1)*max([0 w2bnew(1,1)]) r(1,2)-
2*d(1,2)*max([0 w2bnew(1,2)]) r(1,3)-2*d(1,3)*max([0 w2bnew(1,3)]) r(1,4)-
2*d(1,4)*max([0 w2bnew(1,4)]) -f(1,1)*(1 - e(1,1)) -f(1,2)*(1 - e(1,2)) -f(1,3)*(1 - e(1,3)) -
f(1,4)*(1 - e(1,4)) -f(1,5)*(1 - e(1,5)) -f(1,6)*(1 - e(1,6)) b(1,1) b(1,2) b(1,3) b(1,4) b(1,5)
b(1,6)])*(w2new - w2bnew)' + 0.5* a * (norm(w2new - w2bnew))^2;
clear w2bnew;
Livros Grátis
( http://www.livrosgratis.com.br )
Milhares de Livros para Download:
Baixar livros de Administração
Baixar livros de Agronomia
Baixar livros de Arquitetura
Baixar livros de Artes
Baixar livros de Astronomia
Baixar livros de Biologia Geral
Baixar livros de Ciência da Computação
Baixar livros de Ciência da Informação
Baixar livros de Ciência Política
Baixar livros de Ciências da Saúde
Baixar livros de Comunicação
Baixar livros do Conselho Nacional de Educação - CNE
Baixar livros de Defesa civil
Baixar livros de Direito
Baixar livros de Direitos humanos
Baixar livros de Economia
Baixar livros de Economia Doméstica
Baixar livros de Educação
Baixar livros de Educação - Trânsito
Baixar livros de Educação Física
Baixar livros de Engenharia Aeroespacial
Baixar livros de Farmácia
Baixar livros de Filosofia
Baixar livros de Física
Baixar livros de Geociências
Baixar livros de Geografia
Baixar livros de História
Baixar livros de Línguas
Baixar livros de Literatura
Baixar livros de Literatura de Cordel
Baixar livros de Literatura Infantil
Baixar livros de Matemática
Baixar livros de Medicina
Baixar livros de Medicina Veterinária
Baixar livros de Meio Ambiente
Baixar livros de Meteorologia
Baixar Monografias e TCC
Baixar livros Multidisciplinar
Baixar livros de Música
Baixar livros de Psicologia
Baixar livros de Química
Baixar livros de Saúde Coletiva
Baixar livros de Serviço Social
Baixar livros de Sociologia
Baixar livros de Teologia
Baixar livros de Trabalho
Baixar livros de Turismo
