Planejamento hierárquico da produção em usinas de açúcar e álcool: modelagem de otimização robusta

Hierarchical production planning for sugarcane milling companies: robust optimization model

Planejamento hierárquico da produção em usinas de açúcar e álcool: modelagem de otimização robusta

Hierarchical production planning for sugarcane milling companies: robust optimization model

Rafael Piatti Oiticica de Paiva^I; Reinaldo MorabitoII, ^*

^[e-mail], Usina Santa Clotilde S/A, Brasil

^[e-mail], UFSCar, Brasil

RESUMO

Neste trabalho estuda-se o planejamento hierárquico de produção em usinas de açúcar e álcool e propõe-se um modelo de otimização robusta que considera diversas incertezas nos parâmetros do problema. Esta abordagem pode ser vista como uma alternativa à utilização de programação estocástica robusta para abordar este problema, abordagem esta que foi estudada anteriormente pelos autores. Para resolver os modelos de programação linear e programação inteira mista envolvidos, utiliza-se um software de otimização em programação matemática. Os resultados computacionais obtidos são comparados aos resultados da modelagem determinística de um trabalho anterior dos mesmos autores, utilizando os dados do estudo de caso de uma cooperativa de usinas de açúcar e álcool.

Palavras-chave: Planejamento hierárquico da produção. Setor sucroenergético. Otimização robusta. Planejamento agregado de safra. Análise de incertezas. Programação inteira mista.

ABSTRACT

This work studies the hierarchical production planning of sugarcane milling companies and proposes a robust optimization model that considers several uncertainties in the problem parameters. This approach can be seen as an alternative method for the stochastic robust optimization technique that was previously applied by the authors. Mathematical programming software was used for solving the linear and mixed integer programming problems involved. The computational results obtained are compared with the deterministic approach solutions presented in previous papers using data of a cooperative society of sugarcane milling companies case study.

Keywords: Hierarchical production planning. Sugarcane mills. Robust optimization. Aggregate production planning. Uncertainty analysis. Mixed integer linear programming.

1. Introdução

Recentemente Paiva (2009) e Paiva e Morabito (2012) apresentaram um modelo de otimização determinístico para o problema de planejamento agregado da produção em usinas cooperadas do setor sucroenergético. Esta modelagem, denominada PASUC (Planejamento Agregado de Safra em Usinas Cooperadas do Setor Sucroenergético), considera a relação hierárquica existente entre o planejamento anual da cooperativa e o planejamento tático de safra das usinas cooperadas, propiciando a análise integrada do planejamento de safra em um grupo de usinas. No nível de decisão da cooperativa (primeiro nível), um modelo de programação linear multiprodutos, de dois estágios, multiperíodos, capacitado, com horizonte de planejamento anual e agregação temporal mensal determina as metas de produção de cada usina e define a política de estocagem e o atendimento da demanda. No nível de decisão de cada usina cooperada (segundo nível), um modelo de programação inteira mista multiprodutos, monoestágio, multiprocessos, multiperíodos, capacitado, com horizonte de planejamento do período de moagem (aproximadamente 6 meses) e agregação temporal semanal combina decisões de dimensionamento e sequenciamento dos lotes de produção para o atendimento das metas do primeiro nível hierárquico e define, entre outros, a quantidade de cana-de-açúcar colhida e a quantidade de cana transportada por prestador de serviço e a seleção dos processos de produção de açúcar, álcool, melaço e energia elétrica para cada usina cooperada.

Em artigo posterior, Paiva e Morabito (2011), tendo em vista que o setor sucroenergético possui uma grande quantidade de incertezas que impactam fortemente nesse tipo de planejamento, incorporam incertezas na análise do planejamento de safra utilizando a técnica de programação estocástica robusta (MULVEY et al., 1995; VLADIMIROU; ZENIOS, 1997). Além dessa técnica, outras técnicas de otimização sob incerteza podem ser utilizadas para tratar problemas de planejamento hierárquico da produção, tais como: restrições de chance (CHARNES; COOPER, 1959; JOHNSON; MONTGOMERY, 1974), programação estocástica em dois estágios com recurso (DANTZIG, 1955; BEALE, 1955) e otimização robusta (SOYSTER, 1973; BEN-TAL; NEMIROVSKI, 1998, 1999, 2000; BERTSIMAS; SIM, 2003, 2004; BERTSIMAS; THIELE, 2006). Mais detalhes sobre essas técnicas podem ser encontrados em Joshi (1995), Sox e Muckstadt (1996, 1997), Birge e Louveaux (1997), Sen e Hingle (1999) e Diwekar (2002); além de em pesquisas sobre o estado da arte em otimização sob incerteza (SAHINIDIS, 2004; HERROELEN; LEUS, 2005), programação estocástica (BOUZA, 1993) e otimização robusta (BEYER; SENDHOFF, 2007).

Neste trabalho adota-se a análise de incerteza sob a ótica da otimização robusta proposta por Bertsimas e Sim (2003, 2004) e Bertsimas e Thiele (2006), dado o potencial de aplicação destas técnicas e ao fato de não ser necessário determinar as distribuições de probabilidade que representam as variáveis aleatórias consideradas e, desta forma, facilitar a incorporação da otimização sob incerteza em situações práticas. Propõe-se um modelo de otimização robusta que considera diversas incertezas nos parâmetros do problema e, para resolver os modelos de programação linear e programação inteira mista envolvidos, utiliza-se um software de otimização em programação matemática. Os resultados computacionais obtidos com a abordagem de otimização robusta são comparados aos resultados da modelagem determinística PASUC, utilizando os mesmos dados do estudo de caso da Usina Santa Clotilde (USC) e da Cooperativa Regional dos Produtores de Açúcar e Álcool de Alagoas (CRPAAA), analisados em Paiva e Morabito (2011, 2012). Este artigo está organizado da seguinte maneira: na seção 2 revisa-se resumidamente a abordagem de otimização robusta; na seção 3 aplica-se a abordagem de otimização robusta no problema estudado; na seção 4 analisam-se os resultados desta aplicação no estudo de caso e na seção 5 apresentam-se as considerações finais deste estudo e perspectivas para pesquisa futura.

2. Abordagem de otimização robusta

Alguns autores têm apontado dificuldades para garantir factibilidade e otimalidade na solução de modelos determinísticos e modelos de otimização estocástica em situações onde existe incerteza nos dados de entrada. Exemplos dessas dificuldades são reportados por Ben-Tal e Nemirovski (2000) em testes para avaliar a degradação da solução de modelos de programação linear (PL) quando se impõem pequenas variações (<1%) no valor dos parâmetros de problemas testes disponíveis na NETLIB (www.netlib.org). Uma alternativa para abordar essa questão são os modelos de otimização robusta. Nesses modelos, a garantia de factibilidade é obtida por uma análise de pior caso, que incorpora a incerteza nos parâmetros de forma direta em um modelo de minimização do máximo desvio previsto para as variáveis aleatórias (BERTSIMAS; SIM, 2003). Com essa abordagem de pior caso não é necessário penalizar a infactibilidade do modelo, como em modelos de programação estocástica robusta. Além disso, é possível garantir que a solução obtida é factível para qualquer caso em que os parâmetros incertos variem dentro de um conjunto convexo preestabelecido.

Um dos estudos pioneiros e seminais da teoria de otimização robusta foi apresentado por Soyster (1973), que propôs uma abordagem de pior caso para garantir a factibilidade em modelos de PL com incerteza na matriz tecnológica. Considere o seguinte modelo de PL nas variáveis x_j:

em que cada parâmetro a_ij da matriz tecnológica corresponde ao valor esperado de uma variável aleatória ã_ij, contida em uma hiperesfera com centro em a_ij e raio â_ij , como segue: , m. Note que â_ij corresponde ao desvio máximo da variável aleatória ã_ij com relação ao valor esperado a_ij. Seja η_ij o desvio unitário da variável aleatória ã_ij, definido como: η_ij=(ã_ij - a_ij)⁄â_ij,-[ - 1,1].

De acordo com essas definições, o valor do raio â_ij da hiperesfera determina a inexatidão da variável ã_ij; note que se â_ij= 0 para todo i e j, então se tem um modelo determinístico. Para considerar a incerteza presente nos parâmetros da matriz tecnológica, considere a nova restrição (4) em substituição à restrição (2).

Soyster (1973) propôs uma abordagem determinística de pior caso para garantir a factibilidade da restrição (4) para toda realização â_ij na hiperesfera K_ij. Para isso, considere o seguinte desenvolvimento:

Com essa formulação, o termo determina a "proteção" necessária para garantir a factibilidade de cada restrição i, impondo um espaço de segurança entre e b_i. Substituindo o módulo de x_j por uma variável auxiliar y_j, da seguinte forma: y_j = |x_j|, - y_j< x_j< y_j e y_i≥ 0, com essas novas variáveis e restrições, o modelo (1), (4) e (3) é reescrito por:

Essa formulação agrega robustez ao resultado encontrado, garantindo a factibilidade do modelo para toda realização ã_ij no intervalo K_ij, entretanto também penaliza de forma excessiva a otimalidade da solução (BERTSIMAS; SIM, 2003). Para tratar o conservadorismo do modelo (6)-(9), outros pesquisadores desenvolveram propostas de otimização robusta que agregam maior flexibilidade a esse tipo de estratégia, tais como, Ben-Tal e Nemirovski (1998, 1999, 2000), El-Ghaoui e Lebret (1997), El-Ghaoui, Oustry e Lebret (1998), Bertsimas e Sim (2003, 2004), Bertsimas e Thiele (2006) e Bertsimas, Pachamanova e Sim (2004). Neste artigo, adotam-se as abordagens de Bertsimas e Sim (2003, 2004) e Bertsimas e Thiele (2006), dado a possibilidade de exercer um controle sobre o grau de conservadorismo da solução, entre outras vantagens dessas abordagens.

Considere novamente o modelo de PL definido pelas Equações 1, 4 e 3 com as incertezas do modelo presentes nos coeficientes da matriz tecnológica ã_ij, j - J_i. Por generalidade, J_i é o subconjunto de índices j dos parâmetros da matriz tecnológica que possuem incerteza em cada restrição i, ou seja, não necessariamente todos os parâmetros ã_ij da linha i da matriz tecnológica possuem incerteza. Note que o caso em que todos os parâmetros possuem incerteza corresponde a J_i = {1,2,...,m}. Considere ainda que a distribuição de probabilidade dessas variáveis aleatórias são desconhecidas, porém simétricas nos intervalos [a_ij- â_ij, a_ij + â_ij]. Desenvolvimento análogo a esse pode ser apresentado quando a incerteza está presente no vetor de custos c_j da função objetivo, conforme Bertsimas e Sim (2003). Seja Γ_i o grau de conservadorismo da restrição i, tal que G_i - [0,|J_i|] (por simplicidade, admita que Γ_i seja inteiro) e S_i, um subconjunto de J_i (S_i⊆ J_i) tal que |S_i| = Γ_i. Dessa forma, define-se o seguinte modelo de programação não linear:

em que a função proteção da i-ésima restrição é dada pela seguinte expressão:

Observe nessa formulação que se Γ_i = |J_i| então todas as incertezas da matriz tecnológica são consideradas, o que resulta na formulação de Soyster. Por outro lado, se Γ_i = 0, então nenhuma incerteza é considerada e tem-se um modelo determinístico com os valores esperados a_ij desses coeficientes. Dessa forma, quando se adotam valores inteiros de Γ_i no intervalo [0, |J_i|], tem-se a flexibilidade de ajustar o conservadorismo da solução, como mencionado anteriormente. É interessante notar que a função proteção (14) estabelece um "colchão" de segurança que procura garantir a factibilidade para a i-ésima restrição, independentemente do valor assumido pela variável aleatória que representa a incerteza do problema. Esta garantia de factibilidade cresce à medida que Γ_i se aproxima da cardinalidade do conjunto de parâmetros incertos J_i. Para uma solução factível x_j conhecida, a função de proteção (14) da i-ésima restrição pode ser representada pelo seguinte modelo de otimização (BERTSIMAS; SIM, 2004):

Como esse modelo de otimização é factível e limitado, aplica-se o conceito de dualidade forte e introduzem-se as variáveis duais λ_i e ρ_ij, para que seja possível reescrevê-lo por meio da seguinte formulação:

Para eliminar a função módulo na Equação 19, substitui-se o módulo de x_j por uma variável auxiliar y_j, de forma similar ao que foi feito anteriormente: y_j = |x_j|, - y_j< x_j< y_j e y_j≥ 0. O modelo (18)-(21) é reescrito de forma linear por:

Nesse ponto, é possível substituir a formulação do modelo dual obtido para a i-ésima função proteção b_i(x_j, Γ_i) [Equações 22 a 27] na Equação 12 do modelo não linear original, desconsiderando a minimização por se tratar de um problema de maximização. Com esse procedimento, a formulação obtida conserva a complexidade do modelo original apenas com a adição de algumas variáveis e restrições. A seguir, apresenta-se a formulação do modelo de otimização robusta correspondente:

Uma das vantagens da abordagem proposta por Bertsimas e Sim (2003, 2004) é a possibilidade de inferir, por meio de limitantes, a probabilidade de violação das restrições (4) para os casos em que a variável aleatória ultrapasse o espaço de incerteza considerado. Para calcular o limitante que representa a probabilidade da i-ésima restrição do problema ser violada, considere que x*_j é a solução do modelo de otimização robusta, θ = Γ_i/|J_i|, n = |J_i|, v = (Γ_i + n)/2, e m = v - v. A formulação admite que a distribuição de probabilidade que representa a variável aleatória é simétrica, limitada e independente (BERTSIMAS; SIM, 2004).

Outra opção apresentada por Bertsimas e Sim (2004) para o caso em que , obtida por aproximação de De Moivre-Laplace, pode ser escrita da seguinte forma:

onde Φ é a função distribuição acumulada da normal.

Algumas vantagens da utilização do modelo proposto por Bertsimas e Sim (2004) são: (i) agrega maior flexibilidade por meio do controle do conservadorismo adotado na análise; (ii) não resulta em um modelo de maior complexidade computacional do que o modelo determinístico original (e.g., se o modelo original é um PL, o modelo robusto também vai ser um PL); (iii) pode ser aplicado em problemas práticos com grande quantidade de variáveis reais; (iv) possibilita estimar a probabilidade de violação das restrições; (v) não há necessidade de determinar as distribuições de probabilidade que representam as variáveis aleatórias consideradas. Esse último é aqui considerado como um ponto positivo da metodologia, dado que em geral há dificuldades na definição das distribuições de probabilidade para a situação aqui estudada. Entretanto, outros autores consideram isso como um ponto negativo, por não ser possível incorporar as informações de probabilidade quando as distribuições estão disponíveis (BEYER; SENDHOFF, 2007).

3. Aplicação da abordagem de otimização robusta

Conforme mencionado, os modelos de otimização robusta desenvolvidos e aplicados nesse trabalho baseiam-se no modelo determinístico de dois níveis PASUC, por sua vez composto dos modelos de programação linear PASUC/N1 (nível 1) e de programação inteira mista PASUC/N2 (nível 2). Para que a apresentação dos modelos a seguir fique mais autocontida, reapresentamos no Anexo 1 o modelo PASUC estudado em Paiva (2009) e Paiva e Morabito (2012); para mais detalhes desse modelo, o leitor pode consultar as fontes citadas.

No nível de decisão da cooperativa (PASUC/N1), um modelo de programação linear multiproduto, dois estágios, dinâmico, capacitado, com horizonte de planejamento anual e agregação temporal mensal determina a quantidade e o mix de produção (PU_put) das usinas, a política de estoque (I_pet) e a política de atendimento da demanda (D_pt) da CRPAAA. Essas variáveis devem ser calculadas de forma que a margem de contribuição do grupo seja maximizada e as restrições de mercado, previsão de safra, capacidade de transporte das frotas e capacidade de estoque sejam respeitadas.

No nível de decisão de cada usina cooperada (PASUC/N2), um modelo de programação inteira mista multiproduto, monoestágio, multiprocesso, dinâmico, capacitado, com horizonte de planejamento do período de moagem (aproximadamente 6 meses) e agregação temporal semanal combina decisões de dimensionamento e sequenciamento dos lotes de produção para o atendimento das metas do primeiro nível hierárquico. Esse modelo determina a moagem semanal por tipo de fonte de suprimento (M´_tm)tipo de serviço de transporte f(M"_ft) e por processo k(M""_kt) de forma que a meta de produção do modelo PASUC/N1 seja viabilizada com o mínimo de atraso possível, ao mesmo tempo em que se determina a quantidade de energia elétrica exportada no período t (EE_t). O objetivo é obter um plano tático de safra que proporcione a máxima margem de contribuição, respeitando as restrições de previsão de safra, capacidade de transporte, capacidade de estoque, fluxo de caixa positivo e necessidade de energia térmica e elétrica nos processos industriais.

3.1. Incerteza no preço dos produtos (PASUC/N1-B-VP)

Incertezas em diversos parâmetros do modelo PASUC podem ter impactos importantes na sua solução, entre eles, os valores líquidos de venda de cada produto p em cada período t, i.e., os parâmetros VP_pt da função objetivo do modelo PASUC/N1 (Anexo 1). A incerteza nesse parâmetro é mercadológica em função dos preços do açúcar e álcool estarem sujeitos às variações do mercado de derivativos. Essa característica faz com que seja necessário uma análise complexa sobre os fundamentos do mercado mundial destas commodities, para que seja possível inferir a tendência de variação desse parâmetro. Entretanto, mesmo com a mais completa avaliação não é possível prever com exatidão os preços futuros. Outra característica que torna importante avaliar as variações no parâmetro VP_pt é o forte impacto na função objetivo do modelo PASUC/N1.

O modelo de otimização robusta desta seção é identificado pela abreviação PASUC/N1-B-VP, onde a letra B identifica a utilização da técnica de otimização robusta inicialmente proposta por Bertsimas e Sim (2003, 2004) e as letras VP representam a incorporação da incerteza nos parâmetros VP_pt. Como referência para a modelagem apresentada a seguir, considere as equações do modelo PASUC/N1 (Anexo 1).

Para desenvolver o modelo de otimização robusta com incerteza nos preços dos produtos, define-se a variável aleatória VP~_pt, que representa todos os possíveis valores de preços para o produto p no período t. Considere que o valor médio da variável VP_pt é dado por VP_pt (i.e., o valor nominal do parâmetro no modelo determinístico PASUC/N1) e que o maior desvio não negativo em relação ao valor esperado é dado por VP^_pt (i.e., variação esperado do parâmetro). Dessa forma, tem-se a variável aleatória VP~_pt pertencendo ao seguinte intervalo de incerteza: Seja t^VP o conjunto de coeficientes de todas as variáveis aleatórias VP~_pt que possuem incerteza, G^VP o grau de conservadorismo adotado tal que G^VP - [0,|t^VP|] e S^VP um subconjunto de tal que |S^VP| < G^VP, ou seja, com cardinalidade menor ou igual ao grau de conservadorismo G^VP.

Na sequência, utiliza-se a técnica de modelagem de otimização robusta descrita na seção 2 conforme Bertsimas e Sim (2003, 2004) para reescrever a função objetivo do modelo PASUC/N1. Na formulação que segue, sem perda de generalidade, para todas as variáveis aleatórias do modelo foi considerada apenas a parcela desfavorável dos seus possíveis desvios, ou seja, somente o limite inferior dos seus intervalos de variação . Dessa maneira, obtém-se uma função objetivo alternativa Z^~ no modelo:

s.a. restrições originais PASUC/N1 (Anexo 1)

onde:

e D_pt é uma variável de decisão do modelo que corresponde ao atendimento da demanda (venda) do produto p, no período t (t ou m³) (ver Anexo 1). Suponha agora que a solução ótima da variável D_pt no modelo (39) seja conhecida, digamos D*_pt. Dado que essa variável é não negativa (por definição), pode-se reescrever o segundo termo do lado direito da função objetivo em (39) como o seguinte modelo de PL:

Note que a variável auxiliar z_pt desse modelo não precisa ser inteira, dado que sua solução pode ser obtida por inspeção, simplesmente fazendo-se z_pt= 1 para os G^VP maiores coeficientes na função objetivo. Sejam as variáveis duais λ^VP e associadas, respectivamente, às duas restrições do modelo de PL em (40). Aplicando-se o conceito de dualidade forte, obtém-se a seguinte formulação dual correspondente:

Seguindo o mesmo desenvolvimento apresentado na seção 2, é possível conservar a linearidade do modelo original PASUC/N1 na formulação do modelo de otimização robusta PASUC/N1-B-VP:

s.a. restrições originais PASUC/N1 (Anexo 1)

Note que o modelo PASUC/N1-B-VP tem a função objetivo apresentada em (41) sujeita a todas as restrições do modelo original PASUC/N1 e também às restrições adicionais apresentadas em (42). Em Paiva (2009), também foi estudada outra modelagem de otimização robusta que incorpora a incerteza no parâmetro de eficiência global das usinas ao modelo PASUC/N1 (ver Eatr_ut no Anexo 1). Por motivo de limitação de espaço, esse modelo não será tratado no presente trabalho.

3.2. Incerteza na matriz de rendimentos industriais (PASUC/N2-B-A)

O segundo parâmetro considerado como incerto é a eficiência de conversão de ATR de cada usina, representada pela matriz de rendimentos industriais (A_pkt). As características incertas desse parâmetro são inerentes ao processo produtivo de uma usina de açúcar e álcool, tais como qualidade da matéria-prima processada no período, condições operacionais dos equipamentos industriais, eficiência e eficácia dos operadores e dos insumos utilizados.

Para o desenvolvimento do modelo de otimização robusta com incerteza em A_pkt é necessário reformular o modelo determinístico PASUC/N2 do Anexo 1 para incorporar a variável de produção do produto p no período t (PU_pt), que substitui o termo no modelo original. Essa modificação simplifica a identificação das variáveis de segundo estágio e evita a dependência entre a taxa de produção dos produtos p, em um determinado processo k e período t do parâmetro A_pkt, como pode ser percebido nas equações do modelo PASUC/N2 apresentadas no Anexo 1. Mais detalhes sobre essa correlação estão observados em Paiva (2009).

A dependência faz com que uma das hipóteses adotadas no desenvolvimento da modelagem de otimização robusta em Bertsimas e Sim (2003, 2004) seja violada e impossibilita a aplicação direta dessa técnica ao modelo PASUC/N2 com a formulação original. Entretanto, é possível aplicá-la apenas nos índices de A_pkt que podem ser considerados independentes (k e t) e, dessa forma, incorporar a incerteza na matriz de rendimentos industriais. Essa reformulação pode ser verificada na apresentação do modelo PASUC/N2-B-A adiante, com a substituição do termo do modelo PASUC/N2 do Anexo 1 pela variável PU_pt e a incorporação da restrição (43) descrita abaixo. Diversos testes executados com o modelo PASUC/N2 reformulado comprovaram sua consistência com o modelo PASUC/N2.

Para desenvolver o modelo de otimização robusta PASUC/N2-B-A com incerteza na matriz de rendimentos industriais (A_pkt), define-se a variável aleatória Ã_pkt, onde Ã_pkt - [A_pkt - Â_pkt, A_pkt + Â_pkt]. Seja tª_pt o conjunto de coeficientes k de todas as variáveis aleatórias Ã_pkt que possuem incerteza, Gª_pt, o grau de conservadorismo adotado para cada produto p e período t (neste estudo, adota-se Gª_pt = Gª, ∀p,t), em é um subconjunto de . Seguindo o mesmo procedimento apresentado nas seções 2 e 3.1 para reescrever a restrição (43) do modelo PASUC-N2 (ver Anexo 1), é possível apresentar a seguinte formulação do modelo de otimização robusta PASUC/N2-B-A:

Desta forma, o modelo PASUC/N2-B-A tem a função objetivo do modelo PASUC/N2 modificado (44), incorporando a variável PU_pt para cada subconjunto de produtos p (pa - açúcares, pe - álcoois, PM - melaço), sujeito a todas as restrições do modelo PASUC/N2 [Equações 46-68], além da restrição modificada apresentada em (45), que incorpora as variáveis duais, e as restrições adicionais do modelo de otimização robusta apresentadas em (69).

A formulação proposta para o modelo PASUC/N2-B-A é válida para os índices (processos) k pertencentes ao conjunto tª_pt = {k|Â_pkt > 0}, ∀p,t. Além disso, o número de parâmetros incertos depende do grau de conservadorismo adotado (Gª), podendo variar entre zero (modelo determinístico) e K_pt [modelo de Soyster (1973)], sendo K_pt o número total de processos disponíveis para a produção do produto p no período t.

Entretanto, ao serem analisadas as outras restrições do modelo PASUC/N2-B-A, especificamente as restrições (58) e (59), que tratam da produção tudo ou nada, percebe-se que a variável M""_kt só assume valores diferentes de zero para um único processo k, em cada período t. Sendo assim, o termo que representa a incerteza da matriz de rendimentos industriais no bloco de restrições (69) (Â_pkt M""_kt) é nulo para todos os processos k que não forem selecionados pela variável binária de seleção de processos (x_kt) no período t. Essa particularidade faz com que o modelo PASUC/N2-B-A seja equivalente ao modelo de Soyster (1973) quando Gª for igual a 1, fato que limita a utilização do modelo de otimização robusta acima e torna inócua a utilização de um grau de conservadorismo com valores superiores a 1 (i.e., 1 < Gª < K_pt).

Para contornar essa limitação, pode-se adotar a estratégia proposta em Bertsimas e Thiele (2006) e incorporar o escalonamento da restrição (43), conforme segue:

No texto que se segue, essa segunda abordagem é identificada por PASUC/N2-B-A/A2, enquanto a nomenclatura da primeira abordagem apresentada passa a ser PASUC/N2-B-A/A1. Considere agora que tª_pt representa o conjunto de coeficientes que possuem incerteza em cada uma das (p,t) inequações da restrição (70), que Gª - [0,|tª_pt|], ∀p,t e Sª_pt é um subconjunto de tª_pt (Sª_pt⊆ tª_pt), tal que |Sª_pt| < Gª. A seguir, apresenta-se a formulação do modelo PASUC/N2-B-A/A2 conforme estratégia em Bertsimas e Thiele (2006):

s.a. restrições (46)-(68)

Note no modelo PASUC/N2-B-A/A2 que a restrição (72) incorpora o somatório em t´ ∀p,t. Esse acúmulo de incertezas de forma escalonada propicia o controle do grau de conservadorismo da solução em uma abordagem chamada de quase pior caso (BERTSIMAS; THIELE, 2006), contornando a limitação apresentada pela abordagem de Bertsimas e Sim (2003) para esse tipo de modelo. Em Paiva (2009) também é estudado outro modelo de otimização robusta que considera incerteza no parâmetro de tempo aproveitado de moagem (φ_t) do modelo PASUC/N2. Por motivo de limitação de espaço, esse modelo não é tratado no presente trabalho.

4. Resultados computacionais

Nessa seção são apresentados os resultados computacionais obtidos com a utilização dos modelos de otimização robusta descritos na seção anterior e dos mesmos dados do estudo de caso em Paiva (2009) e Paiva e Morabito (2009, 2012). A seção 4.1 analisa os resultados do modelo PASUC/N1-B-VP e a seção 4.2 analisa os resultados obtidos considerando incerteza tanto no preço dos produtos (VP_pt) quanto na matriz de rendimentos industriais (A_pkt) (PASUC-B). Os resultados aqui reportados referem-se a soluções ótimas e foram obtidos utilizando a linguagem de modelagem GAMS, versão 22.7, com o solver CPLEX 11, em um computador Pentium 4 3,2 GHz, 2GB de memória RAM e sistema operacional Windows XP, service pack 2, ou seja, o mesmo ambiente computacional utilizado por Paiva e Morabito (2012).

4.1. Resultados do modelo de otimização robusta com incerteza nos preços dos produtos (PASUC/N1-B-VP)

Nessa seção são apresentados os resultados obtidos com a utilização do modelo PASUC/N1-B-VP, adotando duas abordagens para tratar a incerteza no preço dos produtos finais (VP_pt). A primeira abordagem (PASUC/N1-B-VP/A1) considera desvios percentuais constantes entre os preços definidos para o modelo determinístico (valor nominal) e o desvio máximo da variável aleatória no modelo de otimização robusta. A segunda abordagem (PASUC/N1-B-VP/A2) considera um crescimento percentual do desvio máximo da variável aleatória. O apelo prático que motivou a utilização dessas duas abordagens foi a suposição de que a incerteza nos preços dos produtos é maior para períodos mais distantes e menor para períodos mais próximos. Dado que o modelo proposto neste trabalho é hierárquico, apresenta-se ao final desta seção uma análise sobre o impacto causado pela incerteza nos preços dos produtos (nível 1) no segundo nível de decisão.

A Figura 1 ilustra a diferença entre as abordagens mencionadas (PASUC/N1-B-VP/A1 e PASUC/N1-B-VP/A2). Os valores dos desvios utilizados para cada uma das abordagens são os mesmos que foram adotados no estudo de caso reportado por Paiva e Morabito (2011) para o cenário A das duas abordagens (PASUC/N1-M-VP/A1 Cen. A, PASUC/N1-M-VP/A2 Cen. A), ou seja, o modelo de otimização robusta aqui descrito utiliza os mesmos desvios percentuais dos cenários pessimistas do modelo de programação estocástica robusta com incerteza no preço dos produtos de Paiva e Morabito (2011).

Figura 1.
Percentual de incerteza em cada abordagem do modelo PASUC/N1-B-VP.

Após a definição das duas abordagens utilizadas para analisar o impacto da incerteza no preço dos produtos, é possível apresentar os resultados obtidos para o modelo PASUC/N1-B-VP, conforme segue. A análise feita nesta seção procura avaliar a deterioração da função objetivo do modelo ao se imporem valores ao parâmetro Γ^VP (parâmetro de controle da robustez do modelo). A variação do parâmetro Γ^VP inicia em zero e segue até o valor em que o aumento de Γ^VP deixa de impactar no resultado do modelo de otimização robusta. A seguir são apresentados os resultados obtidos para a margem de contribuição do modelo PASUC/N1-B-VP (Figuras 2, 3) e do modelo PASUC/N2-B-VP (Figura 4) para cada uma das abordagens.

Figura 2.
Impacto do grau de conservadorismo em PASUC/N1-B-VP/A1 e PASUC/N1-B-VP/A2.

Figura 3.
Impacto do grau de conservadorismo em PASUC/N1-B-VP/A1 e PASUC/N1-B-VP/A2 com G^VP pequeno.

Figura 4.
Impacto do grau de conservadorismo para PASUC/N2-B-VP/A1 e PASUC/N2-B-VP/A2.

A Figura 2 apresenta a curva de deterioração da função objetivo para valores G^VP variando entre 0 e 80 e a Figura 3, os mesmos dados para valores de G^VP variando entre 0 e 10. Por meio dessas duas figuras é possível verificar que a queda percentual da função objetivo é acentuada para valores baixos de G^VP. Ou seja, para uma pequena variação no parâmetro de controle da robustez do modelo, tem-se uma grande variação no valor da função objetivo. Esse resultado indica que existem poucos produtos que representam uma grande parcela da receita das empresas. Por meio dessa constatação é possível dizer que é interessante determinar planos de produção menos dependentes de um único produto ou, caso isso não seja possível, trabalhar com mecanismos de proteção, tais como hedging (mecanismo financeiro que possibilita fixar preços futuros de commodities).

Além disso, percebe-se uma diferença entre as curvas de perda de otimalidade (Figuras 2, 3) para as abordagens PASUC/N1-B-VP/A1 e PASUC/N1-B-VP/A2. Nesse caso, considera-se que a solução apresentada para o modelo robusto da abordagem PASUC/N1-B-VP/A2 antecipa a ocorrência das incertezas dos períodos mais distantes e, consequentemente, propicia uma solução com dependência um pouco menor nos preços desses períodos, fazendo com que o valor de sua solução mais conservadora seja maior que a mesma situação da abordagem A1. O tempo computacional para solucionar os problemas do modelo PASUC/N2-B-VP variou bastante entre as instâncias analisadas (Tabela 1). Todas as soluções obtidas são comprovadamente ótimos globais do problema.

Tabela 1.
Tempo computacional PASUC/N2-B-VP/A1 e PASUC/N2-B-VP/A2 (segundos).

4.2. Resultados do modelo de otimização robusta com incerteza na matriz de rendimentos industriais e no preço dos produtos (PASUC-B)

Nesta seção são apresentados os resultados combinados da incorporação da incerteza no parâmetro de preço do produto (PASUC/N1-B-VP) e também na incerteza proveniente da matriz de rendimentos industriais (PASUC/N2-B-A). No caso da incerteza em A_pkt (nível 2), existem duas abordagens apresentadas (PASUC/N2-B-A/A1 e PASUC/N2-B-A/A2): a primeira está baseada na proposta de Bertsimas e Sim (2003) e a segunda, na modelagem de Bertsimas e Thiele (2006). Conforme discutido anteriormente, o modelo PASUC/N2-B-A/A1 torna-se equivalente ao modelo de Soyster (1973) para valores de Γª maiores que 1 e, portanto, o modelo com maior potencial para incorporar a incerteza na matriz de rendimentos industriais é o PASUC/N2-B-A/A2. Entretanto, testes computacionais executados para esse segundo modelo não obtiveram resultados com gap inferior a 2% em um período de 20 horas de tempo computacional. Levando em consideração as limitações das duas abordagens propostas, optou-se pela apresentação dos resultados de modelo PASUC/N2-B-A/A1.

Antes de iniciar a apresentação dos resultados dessa seção é importante destacar que foram selecionados valores de Γ^VP e Γª de acordo com a probabilidade de violação das restrições do modelo (Tabela 2). Adotou-se Γª=1 para o modelo PASUC/N2-B-A/A1. Vale ressaltar que teste computacionais foram realizados com a utilização de outros valores de Γª e observou-se que não houve alteração nos resultados. No caso de Γ^VP, a Tabela 2 apresenta diversos valores do limitante da probabilidade de violação das restrições de modelo (limitante aproximado; BERTSIMAS; SIM, 2004), em função do conjunto de coeficientes que possuem incerteza (t) e do valor do parâmetro de controle de robustez (Γ), considerando, por exemplo, |t| igual a 7, 24, 144 e 600. Para a aplicação desta seção, adotou-se Γ^VP=20 (|t^VP| = |VP_pt| = 144), de forma que se obtenha uma probabilidade de 5,67% para a violação de cada restrição que possui parâmetros incertos em sua modelagem inicial. A Figura 5 ilustra o comportamento do limitante da probabilidade de violação para diversos valores de t e de Γ.

Tabela 2.
Tabela de limitantes de violação da restrição do modelo de otimização robusta.

Figura 5.
Limitante de violação da restrição Γ^VP e Γª.

Os valores de V^P_pt foram calculados segundo a abordagem 2 do modelo PASUC/N1-B-VP e os valores de Â_pkt, de acordo com um desvio percentual uniforme de 1%. Além disso, utilizou-se a estratégia de redução do número de variáveis binárias existentes no modelo (sem perda de generalidade), conforme discutido em Paiva (2009).

A seguir são apresentados os dados de saída obtidos para o modelo PASUC-B (composição do modelo PASUC/N1-B-VP/A2 e PASUC/N2-B-A/A1), enfocando os resultados de atendimento da demanda (D_pt) por parte da CRPAAA e estabelecimento das metas de produção (Meta_pt) para o caso da USC (primeiro nível); além dos dados de produção (PU_pt), geração de energia (EG_t) e saldo financeiro (S_t) que se referem ao segundo nível.

É interessante comparar os resultados da Tabela 3 com os resultados de Paiva e Morabito (2011, 2012) para avaliar o impacto da incorporação da incerteza na variável de atendimento da demanda (D_pt). Nos resultados apresentados nesses trabalhos, a demanda máxima foi atingida (saldo demanda igual a zero) para o caso dos açúcares Especial, Extra e os alcoóis AEHC e AEAC, diferentemente do resultado apresentado na Tabela 3, onde a demanda máxima de açúcar refinado foi alcançada e a demanda máxima de AEAC não. Esse resultado mostra uma diferença considerável na incorporação da robustez do modelo de programação estocástica robusta (PAIVA; MORABITO, 2011) e a do modelo de otimização robusta apresentado neste trabalho.

Tabela 3. Tabela de atendimento da demanda por produto (t ou m³) - PASUC/N1-B.

A Tabela 4 apresenta a meta de produção mensal da USC (Meta_pt), assim como um comparativo entre a meta acumulada para toda a safra e a produção real da USC na safra 2007/2008. Essa comparação é apenas ilustrativa, pelo fato de ocorrerem alterações normais ao longo da safra no direcionamento do mix de produção, na quantidade de cana processada e em alguns parâmetros considerados nesse processo de planejamento. Os resultados do modelo PASUC/N2-B para a geração de energia estão apresentados na Figura 6, onde se percebe uma geração um pouco superior a 1 milhão Mwh/semana para as primeiras semanas de safra (semana 1 a semana 9) e uma geração máxima de 1.680 Mwh/semana na semana 18.

Tabela 4. Resultado obtido para a meta de produção da USC (t ou m³) - PASUC/N1-B.

Figura 6.
Geração de energia planejada pelo modelo PASUC/N2-B e a dados da safra 2007/2008.

A Figura 7 apresenta a projeção de saldo financeiro ao longo do tempo obtida para o modelo PASUC/N2-B. Analisando essa figura é possível inferir que, nesse cenário, o repasse semanal da cooperativa não está sendo suficiente para cobrir as despesas correntes da usina em análise. Dessa forma, percebe-se um decréscimo do saldo inicial até um valor ligeiramente negativo nas últimas semanas de safra. No caso da última semana de safra (semana 25), verifica-se que o valor do saldo financeiro projetado pelo modelo PASUC/N2 (PAIVA; MORABITO, 2012) é ligeiramente inferior ao projetado pelo modelo PASUC/N2-B (Figura 7). Essa constatação indica que a otimização da margem de contribuição do modelo não está diretamente ligada à maximização do saldo financeiro. Além disso, esse resultado é interessante para embasar negociações com bancos para a obtenção de linhas de crédito de capital de giro e para negociar com a CRPAAA um melhor repasse financeiro (RCa_p, RCb, RCc_t).

Figura 7.
Projeção de saldo financeiro da USC segundo o modelo PASUC/N2-B.

Para finalizar a análise do modelo PASUC-B, são apresentadas duas tabelas de comparação dos resultados globais (Tabelas 5, 6). A Tabela 5 destaca os principais resultados do primeiro nível do modelo PASUC-B (PASUC/N1-B) e faz um comparativo com os resultados obtidos no modelo determinístico (PASUC/N1). Dessa forma, deve-se destacar que os resultados do modelo PASUC/N1-B são idênticos aos resultados do modelo PASUC/N1-B-VP/A2 (seção 3.1) com Γ^VP=20 (Figura 2) e, portanto, as consideração do item 3.1 também são válidas nesse ponto.

Tabela 6. Comparação dos resultados globais - PASUC/N2-B.

Tabela 5.
Comparação dos resultados globais - PASUC/N1-B.

Os resultados da Tabela 6 se referem ao segundo nível de decisão do modelo PASUC-B (PASUC/N2-B). Analisando-se nessa tabela o valor obtido para a margem de contribuição total da USC (função objetivo do modelo PASUC/N2-B), é possível afirmar que o modelo PASUC/N2-B obteve um resultado 12,69% inferior ao valor obtido pelo modelo PASUC/N2, demonstrando o grande impacto decorrente das incertezas no preço dos produtos (nível 1) e na matriz de rendimentos industriais (nível 2) para o planejamento da USC.

Tendo em vista que os modelos de otimização robusta adotam uma abordagem de pior caso, é interessante comparar os resultados da margem de contribuição do modelo PASUC-B (Tabelas 5, 6) com o resultado do cenário mais pessimista (cenário A) do modelo PASUC-M (PAIVA; MORABITO, 2011). Essa comparação indica que o modelo de programação estocástica robusta apresenta valores ligeiramente inferiores aos resultados do modelo de otimização robusta para o caso da margem de contribuição do segundo nível de decisão (R$ 15.833.035 para o cenário A do modelo PASUC/N2-M e R$ 15.915.232 para o modelo PASUC/N2-B). Essa constatação não é observada na comparação da margem de contribuição do primeiro nível (R$ 79.022.260 para os cenários A/D/G do modelo PASUC/N1-M e R$ 85.042.113 para o modelo PASUC/N1-B). Dessa forma, é possível afirmar que o fato do modelo PASUC/N2-B corresponder ao modelo de Soyster (1973) gera um impacto significativo na qualidade da solução do modelo apresentado neste trabalho. O tempo computacional necessário para obter-se a solução ótima do modelo PASUC/N2-B foi 2.314 segundos (~39 minutos).

5. Considerações finais e perspectivas para pesquisa futura

Analisando os resultados obtidos neste estudo é possível verificar que os modelos de otimização robusta permitem analisar o efeito de incertezas existentes nos parâmetros de preço dos produtos, entretanto o mesmo resultado não foi obtido para os parâmetros de rendimento industrial. A seguir são apresentadas algumas considerações sobre a utilização dessa técnica, destacando as vantagens e as limitações percebidas. Algumas vantagens da abordagem de otimização robusta são:

• Não há necessidade de determinar as distribuições de probabilidade que representam as variáveis aleatórias. Quando essas distribuições estão disponíveis, há autores que consideram esse item como ponto negativo, por não ser possível incorporar essas informações de probabilidade (BEYER; SENDHOFF, 2007);

• Possibilita inferir a probabilidade de violação das restrições para os casos em que a variável aleatória ultrapasse o espaço de incerteza considerado;

• Por se tratar de uma abordagem de pior caso, não é necessário utilizar parâmetros de penalização para garantir a factibilidade do modelo;

• Durante a análise da incerteza no preço dos produtos (seção 3.1), possibilita adotar diferentes padrões de comportamento do intervalo de incerteza analisado, tornando flexível a incorporação desse parâmetro do modelo;

• Possibilita exercer controle sobre o grau de conservadorismo da solução, agregando maior flexibilidade de análise e possibilitando a utilização do conceito de diagrama de Pareto robusto (curva de deterioração da função objetivo para diferentes valores Γ). Vale lembrar que esse resultado não foi observado para o caso do modelo PASUC/N2-B-A/A1 e PASUC/N2-B-A/A2, devido a limitações.

Algumas desvantagens da abordagem de otimização robusta são:

• Alta penalização da otimalidade da solução, principalmente quando o número de parâmetros incertos é pequeno. Resultados com perda de otimalidade pequena só foram alcançados com uma alta probabilidade de violação das restrições, ou seja, muito próximo do modelo determinístico (Γ = 0);

• Limitação da aplicação da técnica de otimização robusta em parâmetros que possuem dependência entre seus índices e em restrições de igualdade. Por exemplo, foi necessário uma modificação da modelagem adotada para que a técnica pudesse ser utilizada para a incerteza proveniente da matriz de rendimentos industriais;

• Necessidade da existência de um acúmulo de parâmetros incertos em uma mesma restrição para que o modelo possa controlar a utilização da função proteção e selecionar os parâmetros que mais impactam na função objetivo do modelo. Essa limitação pode ser suplantada pela utilização do modelo de otimização robusta de Bertsimas e Thiele (2006), entretanto, quando o parâmetro incerto é um recurso renovável, a abordagem de Bertsimas e Thiele não pode ser aplicada.

Como sugestão para pesquisas futuras a partir deste estudo destaca-se utilização de modelos de otimização robusta integrados com restrições de chance, como pode ser verificado em Erdoğan e Iyengar (2006) e Chen et al. (2009), ou limitar a análise a um valor de G fracionário entre 0 e 1. Outra possibilidade seria analisar o comportamento dos modelos propostos em termos de desempenho de solução, buscando aceleração da solução dos modelos apresentados por meio de heurísticas baseadas em pacotes de otimização e outros estudos relacionados ao desempenho computacional dessas aplicações. Por fim, sugere-se incorporar incerteza a outros parâmetros dos modelos apresentados neste artigo como, por exemplo, os custos envolvidos, tanto no primeiro nível de decisão quanto no segundo, e melhorar ainda mais a adequação do modelo com a realidade das empresas.

Agradecimentos

Aos revisores anônimos pelos úteis comentários e sugestões, à Usina Santa Clotilde e à Cooperativa Regional dos Produtores de Açúcar e Álcool de Alagoas pelo apoio e pelo fornecimento de dados na realização desta pesquisa. Esta pesquisa contou com apoio do CNPq.

Recebido 05/09/2011;

Aceito 06/05/2012

Anexo 1. Modelos PASUC/N1 e PASUC/N2

Índices do modelo PASUC/N1

u (unidades de produção: usinas cooperadas); t (períodos: meses do ano safra); p (produtos: produtos/coprodutos da fábrica de açúcar e da destilaria de álcool); e (depósitos: conjunto de depósitos próprios e de terceiros)

Parâmetros do modelo PASUC/N1

n_ut Dias disponíveis para moagem em cada usina u e período t (adimensional);

φ_ut Tempo de moagem aproveitado na usina u durante um período t (%);

θ_pa,t Melaço obtido por unidade de açúcar pa produzido na usina u (%);

Matr_u ATR do mel final na usina u (%);

Patr_p ATR correspondente ao produto p (t/t ou m³);

Eatr_ut Eficiência global em ATR na usina u e no período t (%);

ATR_ut ATR da cana da usina u no mês t (kg ATR/t de cana);

Mmin_u Moagem mínima da usina u (t/dia);

Mmax_u Moagem máxima da usina u (t/dia);

Emax_p,e Capacidade máxima de estoque do produto p no depósito e (t ou m³);

CPmax_pu Capacidade de produção diária do produto p na usina u (t ou m³/dia);

CFmax_u Capacidade diária total de produção da fábrica de açúcar na usina u (t/dia);

CDmax_u Capacidade diária total de produção da destilaria de álcool na usina u (m³/dia);

CP_pu Custo de produção do produto p na usina u (R$/t ou m³);

CE_pe Custo de estocagem do produto p no depósito e (R$/t ou m³);

CA Penalização por atraso na entrega da demanda (R$/t ou m³);

Dmin_pt Demanda mínima para o produto p no período t (t ou m³);

Dmax_p Demanda máxima para o produto p durante todo o ano (t ou m³);

VP_pt Valor líquido do produto p no período t (R$/t ou m³);

C_u0 Previsão de moagem de cana para toda a safra da usina u (t);

I_pe0 Estoque inicial de cada produto p no depósito e (t ou m³); EP_p Estoque de passagem de safra do produto p (t ou m³)

Variáveis do modelo PASUC/N1

M_ut Quantidade de cana moída na usina u, no período t (t);

I ⁺_pet Estoque do produto p, no depósito e, no final do período t (t ou m³);

I ^–_pt Atraso na entrega do produto p, no final do período t (t ou m³);

PU_put Produção do produto p, na usina u, no período t (t ou m³);

D_pt Atendimento da demanda (venda) do produto p, no período t (t ou m³)

Conversão de açúcares pa em cristal standard (adimensional);

Função objetivo e restrições do modelo PASUC/N1

A função objetivo acima maximiza a margem de contribuição ao lucro da produção de todos os produtos em todas as usinas e em todos os períodos. As restrições em seguida referem-se, respectivamente: ao balanceamento de estoque de cada produto em cada período; ao atendimento da demanda mínima de cada produto em cada período; à limitação da venda de cada produto na safra ao valor de demanda máxima; à utilização de toda a cana disponível na safra em cada usina, à limitação de moagem de cada usina em cada período; à capacidade máxima de estoque de cada produto em cada depósito e em cada período; à regulação do estoque e atraso de cada produto no último período; às capacidades de produção de cada produto, de todos os produtos da fábrica de açúcar e de todos os produtos da destilaria de etanol, para cada usina em cada período; ao balanceamento de ATR para cada usina em cada período; à utilização do melaço produzido como coproduto da fábrica de açúcar na destilaria anexa de cada usina em cada período; aos domínios das variáveis do modelo.

Índices do modelo PASUC/N2

k (processos dentro da fábrica); t (períodos: semanas de safra que indicam o início e o final da moagem); p (produtos: produtos/coprodutos que podem ser produzidos pela usina); m (matérias-primas: contratos de fornecimento de cana); f (frota: contratos de veículos para transporte de cana)

Parâmetros do modelo PASUC/N2

Conversão de açúcares pa em cristal standard (adimensional);

Conversão dos alcoóis pe em AEAC (adimensional);

Matr_USC ATR do mel final na USC (%); n_t

Dias disponíveis para moagem em cada período t (adimensional);

φ_t Tempo de moagem aproveitado durante um período t (%);

Patr_p ATR correspondente ao produto p (t/t ou m³);

ATR_t ATR da cana no período t (kg ATR/t de cana)

Mmin Moagem mínima da usina (t/dia);

Mmax Moagem máxima da usina (t/dia);

Fmax_f Capacidade máxima de transporte de cana pelos caminhões da frota f (t/dia);

CPmax_p Capacidade de produção diária do produto p (t ou m³/dia);

CFmax Capacidade diária total de produção da fábrica de açúcar (t/dia);

CDmax Capacidade diária total de produção da destilaria de álcool (m³/dia);

α_t Percentual máximo de cana de fornecedores por período t (%);

β_ft Disponibilidade para transporte dos caminhões da frota f, no período t (%);

CF_ft Custo de corte, carregamento e transporte pela frota f, no período t (R$/t);

CA Penalização por atraso na entrega da demanda (R$/t ou m³);

Gfixo Gasto fixo médio em um período t (R$);

GK_kt Gastos variáveis na indústria por processo k, no período t (R$);

GM_mt Gastos variáveis por tipo de matéria-prima m, no período t (R$);

GF_ft Gastos variáveis por tipo de frota f, no período t (R$);

Giro₀

Saldo de caixa inicial da empresa (R$);

RCa_p Adiantamento da CPRAAA pela produção do produto p, no período t (R$/t ou m³);

RCb Adiantamento da CPRAAA pelo ATR equivalente produzido (R$/kg de ATR);

RCc_t Adiantamento extra da CPRAAA obtido no período t (R$/sem)

C₀ Previsão de moagem de cana para toda a safra (t);

Disp_t Parâmetro limitante para o total de cana própria disponível no período t (t);

G Número suficientemente grande (adimensional);

Ib₀ Estoque inicial de bagaço (t);

Fibra_mt Fibra da cana tipo m, no período t (%);

Ub_t Umidade do bagaço após a moenda, no período t (%);

Eb Percentual mínimo de estoque do bagaço produzido (%);

EPb Estoque de bagaço para passagem de safra (t);

RC Rendimento médio das caldeiras (t vapor/t bagaço);

RCF Rendimento médio da casa de força (MWh/t vapor);

CFVAP Consumo fixo de vapor na moagem (t de vapor/t de cana);

CVAP_p Consumo variável de vapor servido em cada produto p (t vapor/t ou m³);

CFE Consumo fixo de energia na moagem (MWh/t de cana);

CVE_p Consumo variável de energia em cada produto p e processo k (MWh/t ou m³);

VAPmax Produção diária máxima de vapor (t/dia);

EGmax Geração diária máxima de energia (MWh/dia);

VE Valor da energia vendida (R$/MWh)

Variáveis do modelo PASUC/N2

X_kt Decisão de utilizar (X_kt = 1) ou não (X_kt = 0) o processo k no período t;

M_t Quantidade de cana moída no período t (t);

Quantidade de cana colhida em cada fornecedor m, no período t (t);

Quantidade de cana transportada pelo tipo de transporte f, no período t (t);

Quantidade de cana moída pelo processo k, no período t (t);

C_mt Quantidade de cana disponível para colheita por matéria-prima m, no período t (t);

S_t Quantidade de capital disponível para giro financeiro no período t (R$);

I ^–_pt Atraso no atendimento da demanda do produto p, no período t (t ou m³);

Ib_t Estoque de bagaço para geração de energia no período t (t);

Mb_t Quantidade de bagaço consumido para geração de vapor no período t (t);

VAP_t Quantidade de vapor produzido no período t (t);

EG_t Quantidade de energia produzida no período t (MWh);

EE_t Quantidade de energia exportada no período t (MWh)

Função objetivo e restrições do modelo PASUC/N2

A função objetivo acima maximiza a margem de contribuição da produção por meio dos processos, da matéria-prima e do tipo de transporte da matéria-prima, em todos os períodos. As restrições em seguida referem-se, respectivamente: ao atendimento da meta de produção acumulada de todos os produtos até certo período; às compatibilidades entre a quantidade de cana colhida, quantidade de cana transportada, quantidade de cana por processo e quantidade de cana moída em todos os períodos; à disponibilidade de cana de cada tipo no início de cada período; à utilização de toda cana disponível para colheita de cada tipo de cana; ao controle de quantidade de cana própria disponível em cada período; ao percentual de cana de fornecedor processada em cada período; à limitação de moagem em cada período; ao balanceamento do saldo financeiro em cada período; à capacidade de transporte de cada frota em cada período; às capacidades de produção de cada produto, de todos os produtos da fábrica de açúcar e de todos os produtos da destilaria de etanol em cada período; à utilização de apenas um processo por período; à determinação em cada período de que a quantidade de cana por processo seja nula sempre que o processo equivalente não esteja sendo utilizado; ao balanceamento de estoque de bagaço em cada período; ao estoque de segurança de bagaço em cada período; ao estoque mínimo de passagem de bagaço no último período de safra; à produção de vapor de acordo com a quantidade de bagaço consumido em cada período; ao balanço de vapor de alta e baixa pressão de toda a planta industrial em cada período; à quantidade de energia excedente que pode ser comercializada em cada período; à capacidade de produção de vapor e energia elétrica em cada período; aos domínios das variáveis de decisão.

• BEALE, E. M. L. On minimizing convex function subject to linear inequalities. Journal of the Royal Statistical Society: Series B, n. 17, p. 173-184, 1955.
• BEN-TAL, A.; NEMIROVSKI, A. Robust convex optimization. Mathematics of Operations Research, v. 23, n. 4, p. 769-805, 1998. http://dx.doi.org/10.1287/moor.23.4.769
• BEN-TAL, A.; NEMIROVSKI, A. Robust solutions of uncertain linear programs. Operations Research Letters, v. 25, n. 1, p. 1-13, 1999. http://dx.doi.org/10.1016/S0167-6377(99)00016-4
• BEN-TAL, A.; NEMIROVSKI, A. Robust solutions of Linear Programming problems contaminated with uncertain data. Mathematical Programming, v. 88, n. 3, p. 411-424, 2000. http://dx.doi.org/10.1007/PL00011380
• BERTSIMAS, D.; PACHAMANOVA, D.; SIM, M. Robust linear optimization under general norms. Operations Research Letters, v. 32, n. 6, p. 510-516, 2004. http://dx.doi.org/10.1016/j.orl.2003.12.007
• BERTSIMAS, D.; SIM, M. Robust discrete optimization and network flows. Mathematical Programming Series B, v. 98, n. 1-3, p. 49-71, 2003. http://dx.doi.org/10.1007/s10107-003-0396-4
• BERTSIMAS, D.; SIM, M. The price of robustness. Operations Research, v. 52, n. 1, p. 35-53, 2004. http://dx.doi.org/10.1287/opre.1030.0065
• BERTSIMAS, D.; THIELE, A. A robust optimization approach to inventory theory. Operations Research, v. 54, n. 1, p. 150-168, 2006. http://dx.doi.org/10.1287/opre.1050.0238
• BEYER, H.; SENDHOFF, B. Robust optimization - A comprehensive survey. Computer Methods Applied to Mechanical Engineering, v. 196, p. 3190-3218, 2007. http://dx.doi.org/10.1016/j.cma.2007.03.003
• BIRGE, J. R.; LOUVEAUX, F. Introduction to Stochastic Programming New York: Springer, 1997.
• BOUZA, C. Stochastic programming: the state of the art. Investigación Operacional, v. 14, n. 2, 1993.
• CHARNES, A.; COOPER, W. W. Chance-constrained programming. Management Science, v. 6, n. 1, 1959. http://dx.doi.org/10.1287/mnsc.6.1.73
• CHEN, W. et al. From CVaR to Uncertainty Set: Implications in Joint Chance Constrained Optimization. Operations Research, v. 58, n. 2, p. 470-485, 2009. http://dx.doi.org/10.1287/opre.1090.0712
• DANTZIG, G. Linear programming under uncertainty. Management Science, v. 1, n. 3, p. 197-206, 1955. http://dx.doi.org/10.1287/mnsc.1.3-4.197
• DIWEKAR, U. Optimization under uncertainty. SIAG/OPT Views-and-News, v. 13, n. 1, p. 1-8, 2002.
• EL-GHAOUI, L.; LEBRET, H. Robust solutions to least-square problems to uncertain data matrices. SIAM Journal on Matrix Analysis and Applications, v. 18, p. 1035-1064, 1997. http://dx.doi.org/10.1137/S0895479896298130
• EL-GHAOUI, L.; OUSTRY, F.; LEBRET, H. Robust solutions to uncertain semidefinite programs. SIAM Journal on Optimization, v. 9, p. 33-52, 1998. http://dx.doi.org/10.1137/S1052623496305717
• ERDOĞAN, E.; IYENGAR, G. Ambiguous chance constrained problems and robust optimization. Mathematical Programming Series B, v. 107, p. 37-61, 2006. http://dx.doi.org/10.1007/s10107-005-0678-0
• HERROELEN, W.; LEUS, R. Project scheduling under uncertainty: Survey and research potentials. European Journal of Operational Research, v. 165, n. 2, p. 289-306, 2005. http://dx.doi.org/10.1016/j.ejor.2004.04.002
• JOHNSON, L.; MONTGOMERY, D. Operations research in production, planning, scheduling and inventory control New York: John Wiley & Sons, 1974.
• JOSHI, R. R. A new approach to stochastic programming problems: Discrete model. European Journal of Operational Research, v. 83, n. 3, p. 514-529, 1995. http://dx.doi.org/10.1016/0377-2217(93)E0218-M
• MULVEY, J. M.; VANDERBEI, R. J.; ZENIOS, S. A. Robust optimization of large-scale systems. Operations Research, v. 43, n. 2, p. 264-281, 1995. http://dx.doi.org/10.1287/opre.43.2.264
• PAIVA, R. P. O. Modelagem do planejamento agregado da produção em usinas cooperadas do setor sucroenergético utilizando programação matemática e otimização robusta 2009. Tese (Doutorado em Engenharia de Produção)-Universidade Federal de São Carlos, São Carlos, 2009.
• PAIVA, R. P. O.; MORABITO, R. An optimization model for the aggregate production planning of a Brazilian sugar and ethanol milling company. Annals of Operations Research, v. 169, p. 117-130, 2009. http://dx.doi.org/10.1007/s10479-008-0428-9
• PAIVA, R. P. O.; MORABITO, R. Programação estocástica robusta aplicada ao planejamento de safra em usinas cooperadas do setor sucroenergético. Gestão & Produção, v.18, n. 4, p. 719-738, 2011. http://dx.doi.org/10.1590/S0104-530X2011000400004
• PAIVA, R. P. O.; MORABITO, R. Otimização do planejamento hierárquico da produção em usinas cooperadas do setor sucroenergético. Produção, 2012. Ahead of print 26 Out 2012. http://dx.doi.org/10.1590/S0103-65132012005000077
• SAHINIDIS, N. V. Optimization under uncertainty: State-of-the-art and opportunities. Computers and Chemical Engineering, v. 28, n. 6-7, p. 971-983, 2004. http://dx.doi.org/10.1016/j.compchemeng.2003.09.017
• SEN, S.; HINGLE, J. L. An introductory tutorial on stochastic linear programming models. Interfaces, v. 29, n. 2, p. 33-61, 1999. http://dx.doi.org/10.1287/inte.29.2.33
• SOYSTER, A. L. Convex programming with set-inclusive constraints and applications to inexact linear programming. Operations Research, v. 21, n. 1, p. 1154-1157, 1973. http://dx.doi.org/10.1287/opre.21.5.1154
• SOX, C. R.; MUCKSTADT, J. A. Multi-item, multi-period production planning with uncertain demand. IIE Transactions, v. 28, n. 11, p. 891-900, 1996.
• SOX, C. R.; MUCKSTADT, J. A. Optimization-based planning for the stochastic lot-scheduling problem. IIE Transactions, v. 29, n. 5, p. 349-357, 1997. http://dx.doi.org/10.1080/07408179708966340
• VLADIMIROU, H.; ZENIOS, S. A. Stochastic linear programs with recourse. European Journal of Operational Research, v. 101, n. 1, p. 177-192, 1997. http://dx.doi.org/10.1016/0377-2217(95)00370-3

*

UFSCar, São Carlos, SP, Brasil

Datas de Publicação

• Publicação nesta coleção
  18 Jun 2013
• Data do Fascículo
  Set 2014

Histórico

• Recebido
  05 Set 2011
• Aceito
  06 Maio 2012