Modelização e Controlo de Sistemas de Manufactura via Filas de Espera

FACULDADE DE ENGENHARIA DA UNIVERSIDADE DO PORTO Modelização e Controlo de Sistemas de Manufactura via Filas de Espera Albino Guambe Mestrado Integrado em Engenharia Eletrotécnica e de Computadores Supervisor: Fernando Lobo Pereira July 29, 2013 c Albino Guambe, 2008 Resumo Esta dissertação endereça um quadro formal de Sistemas com Eventos Discretos para o controlo de sistemas de produção discretos, isto é, constituídos por um número discreto e finito de operações executadas por máquinas, possívelmente distintas, e organizadas de acordo com a especificação do processo de manufactura requerido para obter os produtos acabados inerentes à especificação do sistema. Este quadro formal permite a modelização do processo de produção discreto como uma rede de filas de espera. Para fundamentar este facto, foi efectuada uma revisão das bases fundamentais subjacente à modelização desta classe de sistemas focadas na teoria dos autómatos de estado finitos, temporizados e estocásticos e, mais precisamente, nos sistemas com a propriedade de Markov. A razão da escolha desta classe assenta nos factos de ser computacionalmente tratável e da gama de aplicações práticas endereçável neste contexto ser extremamente vasta. Efectuamos uma revisão dos principais conceitos e resultados relativos a filas de espera, tendo esta dissertação focado na classe M/M/1 que é a mais simples. A razão para esta opção consistiu em tornar clara a perspectiva sistémica da metodologia geral inerente a este quadro formal e, ao mesmo tempo, ilustrar de forma convincente a efectividade desta através do tratamento completo de um exemplo simples apresentado no penúltimo capítulo desta dissertação. Finalmente, importa observar que não só se procurou deixar clara a potencialidade da generalidade das bases fundamentais que permitem endereçar classes de sistemas mais genéricas do que a considerada aqui, como também, para este objectivo, se apontaram alguns desafios e caminhos para os abraçar. i ii Abstract This dissertation concerns a formal framework for Discrete Event Systems for the control of discrete production systems, that is, systems composed by a discrete, finite set of operations executed by, possibly distinct, machines, and organized according the specification of the manufacturing process required to obtain the finished products inherent to the specification of the system. This formal framework enables the modeling of the discrete production system as a network of queues. To show this fact, a review of the fundamental foundations underlying the modeling of this class o systems is provided with a focus on automata theory, encompassing finite state, timed and stochastic automata, and, more precisely, on the systems exhibiting the Markov property. The reason for the choice of this class of systems relies on the facts that it is highly computationally tractable and that its range of practical applications is extremely wide. We reviewed the key concepts and results concerning queueing theory, with a special focus on the M/M/1 class, which is the simplest. The key reasons for this option consisted in clarifying the systemic perspective of the general methodology inherent to the adopted formal framework, and, at the same time, illustrating in a convincing way its effectiveness by fully treating a simple example that we present in the chapter before the last one. Finally, it is important to observe that, we sought not only to let it clear the generalization potential of the framework foundations that allow to treat classes of systems much wider than the one presented here, but also to point out some challenges and research avenues to pursue them inherent to this desideratum. iii iv Agradecimentos Esta dissertação é o culminar de um longo percurso. Tão longo que se torna difícil encontrar palavras de agradecimento para com a Instituição, a FEUP. Não posso deixar de incluir uma palavra de sincero agradecimento aos muitos amigos - colegas, docentes e funcionários - que me ajudaram a tornar este percurso mais curto. Sendo tantos, não vou tentar nomeá-los a todos sob pena do esquecimento me levar a incorrer nalguma injustiça profunda. Contudo, vou abrir uma excepção, e deixo aqui uma palavra de profunda gratidão ao meu amigo Zé Bento sem o qual possivelmente esta dissertação não estaria a ser defendida hoje e aqui. Por fim, deixo a palavra de agradecimento mais importante para o futuro, a minha família, que, com infinita paciência e amor, sempre me apoiou neste longo e difícil percurso. Porto, 25 de Junho de 2013 Albino Guambe v vi Contents 1 Introdução 1 2 Revisão sobre Autómatos 5 2.1 AutómatosdeEstado................................ 5 2.1.1 Linguagens................................. 7 2.1.2 Operações sobre Autómatos . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2 Autómatos Temporizados e Autómatos Estocásticos . . . . . . . . . . . . . . . . 12 2.2.1 Autómatos Temporizados . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2.2 Autómatos Temporizados e Estocásticos . . . . . . . . . . . . . . . . . . 14 3 Processos de Markov 19 3.1 ProcessosdePoisson ................................ 19 3.2 CadeiasdeMarkov ................................. 22 3.3 Resposta de Cadeias de Markov . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.4 Cadeias de Markov em Tempo Contínuo . . . . . . . . . . . . . . . . . . . . . . 28 4 Teoria das Filas de Espera 33 4.1 Especificação de Filas de Espera . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.2 Medidas de Desempenho e Dinâmica de Filas de Espera . . . . . . . . . . . . . 35 4.3 FiladeesperaM/M/1................................ 37 5 Aplicação ao Controlo e Optimização de Sistemas de Manufactura Discretos 41 5.1 Modelização e Análise de um Sistema de Manufactura . . . . . . . . . . . . . . 41 5.2 Redes de Filas de Espera Markovianas . . . . . . . . . . . . . . . . . . . . . . . 43 5.3 Resoluçãodoexemplo ............................... 46 6 Conclusões e Trabalho Futuro 49 References 51 vii 4Introdução Chapter 2 Revisão sobre Autómatos Neste capítulo apresentamos um estado da arte seleccionado sobre modelização e controlo de autómatos, uma vez que constitui o paradigma escolhido para modelizar o comportamento de sistemas de manufactura, a classe de aplicações considerada nesta tese que, como vimos no capítulo 1, é convenientemente modelizada como Sistemas com Eventos Discretos (SED). Seguimos de perto a referência [1]. Deste modo, na próxima secção, apresentamos uma revisão sobre autómatos finitos, incluindo algumas noções de linguagens regulares que lhe estão associadas, bem como operações com autómatos que são importantes para modelização, análise e síntese desta classe de sistemas. Em particular, estas operações são relevantes para a síntese de sistemas de controlo supervisor. Em determinadas circunstâncias, é possível controlar o comportamento - quer em regime transitório quer em regime permanente - de um SED modelizado por autómato através de operações com outros autómatos convenientemente projectados como controladores, [6,3]. O âmbito desta dissertação não permite abordar com um mínimo de propriedade as questões extremamente importantes que se colocam em controlo supervisor. Assim, focaremos apenas na determinação de regimes estacionários para sistemas modelizados por certas classes de autómatos temporizados e estocásticos, uma vez que as operações dos sistemas de manufactura demoram tempo e estão sujeitos a perturbações de carácter aleatório. Estas são as razões para incluir estas extensões da Teoria dos Autómatos numa última secção. De forma sucinta, consideraremos primeiro autómatos temporizados e, seguidamente, estenderemos estes para o contexto estocástico que surge não só na variável tempo associada à ocorrência de eventos como também no mecanismo de transição de estado. 2.1 Autómatos de Estado Um Sistema com Eventos Discretos (SED) é um sistema com um número de estados discretos, cuja evolução do estado depende da ocorrência de eventos (“event-driven”) assíncronos que ocorrem ao longo do tempo. Assim, a um evento corresponde sempre uma transição de estado - nem que seja para o mesmo - bem como, eventualmente, uma mudança de variáveis associadas. No que respeita às transições de estado num SED, as principais suposições são: 5 6Revisão sobre Autómatos •qualquer transição é despoletada instantaneamente, •ocorre apenas uma transição num dado instante, e •o número de transições é finito. Como referimos na introdução, o conceito de SED é muito popular na modelização de sistemas que surgem em processos de manufactura, sistemas computacionais, de comunicação, rodoviários, ou seja, sistemas que envolvem de alguma forma filas de espera, bem como em qualquer classe de sistemas referidos como híbridos. Um abordagem muito comum para representar o comportamento de SED’s consiste na utilização de autómatos, determinísticos ou não determinísticos. Um autómato determinístico consiste num sextúplo definido por, [1,7]: G= (X,E,f,Γ,x0,Xm),(2.1) onde 1. Xé o espaço dos estados discreto, 2. Eé um conjunto de eventos finito, 3. f:X×E→Xé a função de transição de estado, 4. Γ:X→2Eé o mapa de activação de eventos em cada estado, ou seja, Γ(x) = {e∈E:f(x,e)está bem definido}, 5. x0é o estado inicial, e 6. Xm,Xm⊂Xé o subconjunto de estados marcados, ou seja, os estados para os quais se considera que o sistema cumpriu um dos seus objectivos. Qualquer autómato pode ser representado graficamente por um diagrama de transição de estados, e, consequentemente, o comportamento de um SED pode ser especificado graficamente usando o diagrama de transição de estados. Um exemplo trivial de um SED é dado na figura 2.1. A figura 2.2 apresenta um autómato que exibe não-determinismo. A ocorrência do evento a pode dar origem a suas transições distintas sem que qualquer informação permita determinar qual é que pode ocorrer. Tal como em sistemas físicos tradicionais, por exemplo, cujo estado evolui integrando uma equação diferencial, é possível considerar autómatos com entradas e saídas. Existem duas possibilidades para o efeito: •Autómato de Moore. A cada estado corresponde uma saída que pode ser vista como uma extensão da marcação. Neste caso, o autómato pode ser considerado como tendo duas saídas: uma marcada e outra não marcada. 2.1 Autómatos de Estado 7 Figure 2.1: Exemplo de um autómato simples •Autómato de Mealey. Cada transição é etiquetada por um evento da forma “evento de entrada/evento de saída”. Tais eventos especificam quais as entradas que podem ocorrer num dado estado e qual a saída que o autómato produz quando muda de estado. A figura 2.3 ajuda a ilustrar as diferenças entre estes dois tipos de autómatos, através da representação de um dado sistema em cada um deles. Uma das vantagens desta representação de SED é que estão disponíveis muitas técnicas de Teoria dos Automátos e das linguagens para apoiar a análise destes sistemas. 2.1.1 Linguagens Nesta subsecção, apresenta-se uma revisão de alguns conceitos de linguagens (ver [8] para mais detalhes). O conjunto de eventos Ede um SED designa-se de alfabeto, onde Eé um conjunto finito de símbolos, sendo uma sequência de símbolos sde um alfabeto designada de “string”. |s|denota o comprimento da “string” se o símbolo εrepresenta uma “string” vazia. Uma linguagem Ldefinida em Eé um conjunto de “strings” de Ecom comprimento finito. Seja E={a,b,g}, então exemplos de linguagens em E, são •L1={ε,a,abb}, e •L2={todas as “strings” possíveis de comprimento 3 iniciadas com o evento a}. Em termos práticos, uma linguagem corresponde a um conjunto de sequências de eventos admissíveis que correspondem a dados conjuntos de operações a ter lugar num sistema real representado no contexto formal de SEDs. Para a análise destes sistemas importa pois definir um conjunto de conceitos e de operações sobre linguagens. O fecho de Kleene (Kleene-closure) denota-se por E∗e define-se como o conjunto de todas as “strings” finitas mas de comprimento arbitrário de símbolos de E. Por exemplo, se E={a,b,c}, então E∗={ε,a,b,c,aa,ab,ac,ba,bb,bc,ca,cb,cc,aaa,...}. 8Revisão sobre Autómatos Figure 2.2: Exemplo de um autómato exibindo não determinísmo As noções de sub-“string”, prefixo e sufixo são uma parte importante da terminologia. Seja tuv =s, então té prefixo de s,uis sub-“string” de sevis sufixo of s. Importa definir algumas operações sobre linguagens, [8,1]. Considere-se E={a,b,g},L1= {ε,a,abb}eL4={g}nos exemplos abaixo para ilustrar as operações. •Concatenação: Sejam La,Lb⊆E∗, então LaLb={s∈E∗:s=sasb,sa∈Laesb∈Lb}. Exemplo: L1L4={ε,ag,abbg}. •Fecho-em relação aos prefixos: Seja L⊆E∗, então ¯ L:={s∈E∗:∃t∈E∗tal que st ∈L}, ou seja, ¯ Lé o conjunto de todos os prefixos de L, e L⊆¯ L. Exemplos: ¯ L1={ε,a,ab,abb}e¯ L4={ε,g}. •O fecho de Kleen denota-se de L∗. Seja L⊆E∗, então L∗:={ε}[L[LL[.... Exemplos: L∗ 1={ε,a,abb,aa,aabb,...}eL∗ 4={ε,g,gg,ggg,...}. •Pós-linguagem. Dados L⊂E∗es∈L,L/sdesigna a Pós-linguagem de Ldepois de s e define-se por L/s:={t∈E∗:st ∈L}. Obviamente, L/s=/0 sempre que s/∈¯ L. Exemplos: L1/ab ={bg}eL4/g={ε}. É importante observar que ε6=/0 e que {ε}é uma linguagem não vazia uma vez que contém a “string” vazia. Por outro lado, se L=/0, então ¯ L=/0, e se L6=/0, então ε∈¯ L. Finalmente, não é difícil concluir que /0∗={ε}e que {ε}∗={ε}. Os operadores projecção e projecção inversa são extremamente importantes na análise de sistemas com eventos discretos. Sejam E1eE2dois conjuntos de eventos com E2⊂E1. O operador 2.1 Autómatos de Estado 9 Figure 2.3: Exemplo de um autómato com entrada e saída nas formas de Moore e de Mealey projecção P:E∗ 1→E∗ 2define-se, naturalmente, como sendo P(ε) = ε P(e) = (ese e∈E2 εse e∈E1\E2 P(se) = P(s)P(e)para s∈E∗ 1,e∈E2. Naturalmente, a projecção inversa é um operador que não está unicamente definido para cada “string” e, portanto toma como valor conjuntos. Tem-se que P−1:E∗ 2→2E∗ 1é dado por P−1(t) = {s∈E∗ 1:P(s) = t}. Uma linguagem pode ser considerada como um modo formal de descrever um SED, uma vez que especifica todas as sequências admissíveis de eventos que um SED pode gerar. Por outro lado, um autómato é um dispositivo capaz de representar uma linguagem através de regras bem definidas especificada pelo diagrama de transição de estado. A linguagem gerada pelo autómato G= (X,E,f,Γ,x0,Xm)é dada por L(G) = {s∈E∗:f(x0,s)está defindo}. A linguagem marcada pelo autómato Gé definida por Lm(G) = {s∈L(G):f(x0,s)∈Xm}. Sendo os autómatos dispositivos práticos muito convenientes para simular e analisar SEDs, uma questão importante que se coloca é se é sempre possível representar uma dada linguagem por um autómato. Observe-se que a linguagem poderá em muitas situações requerer “strings” com um número infinito de símbolos. A resposta a esta questão é positiva, embora possa requerer autómatos com um número infinito de nós. Contudo, importa observar que é possível gerar linguagens incluindo “strings” com um número infinito de símbolos com autómatos de estado finitos. As linguagens que podem ser geradas por autómatos de estado finitos dizem-se regulares. Prova-se que a classe de linguagens geradas por autómatos de estado finitos exibindo determinismo é a mesma que a dos autómatos de estado finito exibindo não determinísmo. 10 Revisão sobre Autómatos Prova-se o seguinte resultado importante: Sejam L1eL2linguagens regulares, então ¯ L1,L∗ 1,Lc 1(= E∗\L1),L1[L2,L1L2,eL1\L2 também são regulares. 2.1.2 Operações sobre Autómatos A motivação para considerar operações sobre autómatos reside no facto de serem extremamente relevantes nas seguintes actividades de projecto de SEDs: •Modelização - Representar autómatos complexos combinando ou compondo outros autómatos mais simples, bem como obter representações mínimas, simplificando os autómatos através da eliminação dos estados não relevantes. •Análise do comportamento de autómatos - Determinar se a resposta do sistema possui as propriedades desejadas, nomeadamente, admissibilidade dos objectivos definidos, estabilidade, desempenho, segurança, e outros requisitos pré-especificados para o sistema. •Projecto de Sistemas - Modificar autómatos, tendo em vista obter uma resposta do sistema (quer em regime transitório quer em regime permanente) que correspondam a comportamentos pré-especificados ou a optimizar um dado critério de desempenho do sistema. Por exemplo, controlar um autómato pode ser obtido efectuando a composição do autómato dado representando o sistema original com outro designado de controlador. As operações sobre autómatos podem ser classificadas de unitárias, ou seja que envolvem apenas um autómato, ou actuarem em mais de que um autómato. Consideraremos apenas o conjunto de operações mais significativas e, estabeleceremos as relações com as linguagens associadas. Para mais detalhes, considere-se [1,5,7]. As operações ditas unárias mais importantes são as seguintes: •Ac(A)- Extracção da Parte Acessível do autómato A. Esta operação consiste em eliminar os estados que não são necessários e as transições associadas. Os estados que não são necessários são aqueles que nunca são atingidos. Ac(A)também é designada da parte ou componente atingível de A. Observe que esta operação não tem qualquer efeito sobre a linguagem L(A)(ou a linguagem marcada Lm(A)) associada ao autómato. Esta operação é importante para simplificar o autómato sobretudo após operações de composição. •CoAc(A)- Co-Acessibilidade. Eliminação de todos os estados do autómato Apara os quais não existe uma sequência de eventos que levem a qualquer um dos estados marcados, ou seja, em Xm. Um estado qde um autómato Aé dito Co-Acessível se existe um sequência de eventos sque levam o sistema de qaXm, ou seja, f(q,s)∈Xm. Esta operação pode afectar 2.1 Autómatos de Estado 11 L(A)mas não altera Lm(A). Se A=CoAc(A)então diz-se que o autómato é Co-Acessível. Observe que se um autómato é não bloqueante, então tem de ser Co-Acessível. •Trim(A)- Operação “Trim”. Esta operação define-se a partir das anteriores do seguinte modo: Trim(A):=CoAc(Ac(A)) = Ac(CoAc(A)), sendo a ordem irrelevante. O autómato Adiz-se “Trim” se for simultaneamente Acessível e Co-Acessível. •Ac=Comp(A)- Complemento do autómato A. Considere que Aé um autómato “Trim”, então Acé um autómato que “marca” a linguagem E∗\L(A)definido por X[{xd},E,ft,x0,X[{xd}\Xm, sendo xdum estado não marcado para “dumping”, e fto completamento de fpor forma a ter f(x,e) = xdpara todas as transições f(x,e)que não estejam definidas e, em particular, f(xd,e) = xdpara qualquer evento e∈E. Note-se que, como seria de esperar, se tem L(Ac) = E∗eLm(Ac) = E∗\Lm(A). As operações envolvendo vários autómatos que consideraremos são as que correspondem aos dois tipos de composição de dois autómatos. Estas operações podem ser estendidas de forma óbvia a múltiplos autómatos. A relevância destas operações prendem-se com o facto de que o projecto de um controlador para um SED representado por um autómato consiste em construir um autómato controlador que, uma vez composto com o autómato dado, produz um novo autómato com as propriedades (resposta - quer nos regimes transitório quer permanente - e desempenho) pretendidos para o sistema controlado. O projecto de controladores designa-se na literatura de controlo supervisor, [3,6]. Considerem-se os autómatos AeB, dados, respectivamente, por A= (XA,EA,fA,ΓA,xA,0,XA,m) eB= (XB,EB,fB,ΓB,xB,0,XB,m), onde, para facilitar, consideramos ΓA(x) = EA,∀x∈XA, e ΓB(x) = EB,∀x∈XB. •Composição Produto: A×B:=Ac(XA×XB,EATEB,f,Γ,xA,0·xB,0,XA,m×XB,m)onde a operação ×entre dois conjuntos denota o seu produto Cartesiano, a operação ·de dois estados denota a sua ordenação num par, e –f(xA·xB,e) = fA(xA,e)·fB(xB,e)se e∈ΓA(xA)∩ΓB(xB)e está indefinida no caso de e/∈ΓA(xA)∩ΓB(xB). –XA,m×XB,mé o produto cartesiano de todos os estados marcados. Um par de estados em que um deles não é marcado, é considerado não marcado. Repare-se que eventos são considerados sincronizados para ambos autómatos e ocorrem sempre que ocorrem para ambos. Decorre daqui que L(A×B) = L(A)∩L(B)e que Lm(A×B) = Lm(A)∩Lm(B). 12 Revisão sobre Autómatos •Composição Paralela: AkB:=Ac(XA×XB,EASEB,f,Γ,xA,0·xB,0,XA,m×XB,m). As operações ×entre conjuntos e ·entre estados é como acima. Agora, tem-se f(xA·xB,e):=           fA(xA,e)·fB(xB,e)se e∈ΓA(xA)TΓB(xB) fA(xA,e)·xBse e∈ΓA(xA),e/∈ΓB(xB) xA·fB(xB,e)se e/∈ΓA(xA),e∈ΓB(xB) não definido se e/∈ΓA(xA),e/∈ΓB(xB) Observe que se EA=EBeΓA=ΓB, então esta composição coincide com a composição produto. 2.2 Autómatos Temporizados e Autómatos Estocásticos É do conhecimento geral que muitos sistemas - como por exemplo, sistemas de transporte, sistemas de comunicação, etc. -, que operam em ambientes que estão sujeitos a incertezas de forma continuada. Por exemplo, uma simples máquina envolvendo apenas um botão a ser operado por seres humanos pode ter um comportamento imprevisível pelo facto do comportamento humano ser imprevisível ou, simplesmente, devido à ocorrência de avarias imprevisíveis. 2.2.1 Autómatos Temporizados A variável tempo desempenha um papel muito importante no controlo de sistemas de manufactura: as operações requerem tempo para execução e o custo das operações de manufactura dependem fortemente da duração do período de tempo que os produtos levam a ser manufacturados. Este incluí não só a duração das operações em si mas também os tempos de espera antes destas (inventário distribuído). No caso de SED, o tempo é introduzido através da temporização dos eventos, [9]. Ou seja, para cada tipo de eventos, define-se um relógio que determina o próximo instante em que o evento vai ocorrer. Nesta subsecção, consideramos que o relógio está completamente pré-determinado e que os seus elementos são determinísticos e seguiremos de perto a notação em [1]. Seja E={ei:i=1,...,NE}um conjunto de eventos finitos. Um “relógio”, ou estrutura temporal, é um vector de sequências, possivelmente infinitas, de números positivos que, para cada tipo de evento, indica a duração do intervalo entre duas ocorrências sucessivas desse evento, ou seja, V={{vi,1,vi,2,...,vi,k,...}:ei∈E,vi,k>0}. Se ti,k,k=0,1,2,..., denotar os instantes de ocorrência do evento ei, então, vi,k=ti,k−ti,k−1, k=1,2,.... As seguintes noções são importantes para explicar como é que o relógio determina o próximo evento a decorrer. Observe-se que o processo fica mais complicado pelo facto de, dependendo do estado activo do autómato, um dado evento poder estar inibido de ocorrer. Seja ttal que ti,k−1<t<ti,k, então 2.2 Autómatos Temporizados e Autómatos Estocásticos 13 •yi,k=ti,k−té o tempo de vida residual do evento eina k-ésima ocorrência de eventos, ou o valor corrente do relógio para o evento ei. •zi,k=t−ti,k−1é a idade do evento eina k-ésima ocorrência de eventos. Observe que, obviamente, se têm vi,k=yi,k+zi,k. •Ni,kdesigna o número de ocorrências do evento ei∈Eno intervalo de tempo [t0,tk]onde tk é o instante da k-ésima ocorrência de eventos no sistema correspondende a ktransições de estado. Para descrever o funcionamento do relógio necessitamos da seguinte notação: •(x,e,t)é um triplo formado pelo estado actual, x, o último evento ocorrido, e, (ou seja, o que levou ao estado actual) e to instante mais recente em que o evento eocorreu. •Para cada evento ei, considera-se Nieyicomo sendo, respectivamente, o número de ocorrências do evento eiaté ao momento actual, e o valor corrente do relógio para o evento ei. •(x0,e0,t0)é um triplo formado pelos estado, evento e instante de ocorrência seguintes. Obviamente, ter-se-á de ter x0=f(x,e0),e0∈Γ(x)et0é o instante em que e0vai ocorrer. •De forma semelhante, também se definem N0 iey0 i. O mecanismo do relógio consiste pois em determinar o triplo (x0,e0,t0)que vai constituir a próxima transição de estado, bem como calcular os valores de N0 iey0 iassociados. Uma implementação natural do relógio pode ser efectuada pelo seguinte algoritmo, [1], 1 Dado o estado x, determinar o conjunto de eventos admissíveis Γ(x). 2 Para cada evento ei∈Γ(x), calcular yie determinar o evento a que corresponde o menor valor do relógio y∗, ou seja, y∗=min ei∈Γ(x){yi}. 3 O próximo evento e0a ocorrer é aquele a que corresponde o valor de relógio y∗, ou seja, e0=arg min ei∈Γ(x){yi}. 4 Calcular o novo estado: x0=f(x,e0). 5 Calcular o instante da próxima ocorrência do evento e0:t0=t+y∗. 6 Actualizar os novos valores do relógio para os eventos admissíveis em x0, ou seja, para ei∈Γ(x0): y0 i=(yi−y∗se ei6=e0eei∈Γ(x) vi,Ni+1se ei=e0ou ei/∈Γ(x). 20 Processos de Markov b) As variáveis aleatórias N(t),N(t,t1),N(t1,t2),...,N(tk−1,tk), são mutuamente independentes para 0 ≤t≤t1≤... ≤tk. Esta hipótese declara que os incrementos do processo são independentes, ou seja, que não existem “interferências” entre os diferentes intervalos de tempo no que respeita a ocorrências de eventos. c) A probabilidade P[N(tk−1,tk) = n],n=0,1,..., depende de tk−tk−1mas não de tk−1etk em si. Este axioma, juntamente com o anterior, afirma que o processo de contagem tem incrementos independentes e estacionários. Observe-se que, mesmo que alguns destes axiomas não sejam satisfeitos, por exemplo, existência de períodos preferênciais, ocorrência mais do que um evento simultaneamente, etc., é, muitas vezes, possível redefinir o sistema por forma a que os referidos axiomas sejam satisfeitos. Destes axiomas, demonstra-se que Pn(t)satisfaz a distribuição de Poisson, ou seja, existe um número λ>0 tal que, para n=0,1,2,... e para t>0, se tem Pn(t) = (λt)n n!e−λt. O valor esperado e a variância podem ser calculados pelas correspondentes fórmulas gerais: •Valor esperado E[N(t)] = ∞ ∑ n=0 nPn(t) = e−λt(λt) ∞ ∑ n=1 (λt)n−1 (n−1)!=λt. •Variância Var[N(t)] = E[(N(t)−λt)2] = E[N(t)]2−(λt)2=λt, uma vez que E[N(t)]2= ∞ ∑ n=0 n2Pn(t) = e−λt ∞ ∑ n=1 n(λt)n−1 (n−1)!= [(λt)2+(λt)]e−λt ∞ ∑ n=0 (λt)n n!. Dada a importância dos processos de Poisson no estudo das Cadeias de Markov, vamos estudar as suas propriedades. Começamos por considerar o caso em que os eventos que ocorrem são de um único tipo. •Determinação da distribuição da variável duração do intervalo de tempo entre eventos (ou mais simplesmente “tempo entre eventos”), Gk(t) = P[Vk≤t]. Assuma-se que o evento k−1 ocorreu em Tk−1=tk−1. Como P[Vk>t] = 1−Gk(t)e P[Vk>t|Tk−1=tk−1] = P[N(tk−1,tk−1+t) = 0] = P0(t) = e−λttem-se que Gk(t) = 1−e−λt. Daqui se conclui rapidamente a propriedade de ausência de memória, uma vez que: P[V≤z+t|V>z] = P[z≤V≤z+t] P[V>z]=(1−e−λ(z+t))−(1−e−λz) e−λz=1−e−λt=P[V≤t]. 3.1 Processos de Poisson 21 •Uma variável importante é o tempo de vida residual, designada de Y, ou seja, o tempo que falta para a próxima ocorrência, interessando pois determinar a sua função de distribuição de probabilidade. Para uma dada estrututra de relógio V, o tempo de vida residual é, naturalmente, referida à idade do evento que se verifica no instante actual, z. Seja H(t,z) = P[Y≤t|V>z]. Logo, uma vez que V=Y+z, H(t,z) = P[V≤z+t|V>z] = 1−e−λt. Importa observar o facto da função distribuição de probabilidade do tempo de vida residual não depender da idade do evento. Assim, para o processo de Poisson de parâmetro λ, está automaticamente associada uma distribuição exponencial para a duração dos intervalos entre-eventos, G(t) = 1−e−λt, e, como consequência, resulta a propriedade de ausência de memória, ou seja, o facto da função distribuição de probabilidade do tempo de vida residual H(t,z) = G(t)não depender da idade z do evento. Importa agora, verificar como é que as propriedades anteriores se estendem para o caso em que m,m>1, tipos de eventos poderão ocorrer, ou seja, são viáveis no estado corrente do sistema, que, em conformidade com os axiomas acima, são processos mutuamente independentes. Já verificamos que o próximo tempo “entre eventos” é o menor dos tempos de vida residuais, ou seja Y∗=min i=1,...,m{Yi}onde cada variável aleatóriaYitem uma função distribuição de probabilidade Gi(t) = 1−e−λit. Assim, usando o facto de que [min i=1,...,m{Yi}>t] = ∩m i=1[Yi>t], bem como a independência dos processos, temos que P[Y∗≤t] = 1−P[Y∗>t] = 1−P[min i=1,...,m{Yi}>t] = 1− m ∏ i=1 P[Yi>t]. Fazendo Λ= m ∑ i=1 λi, temos que P[Y∗≤t] = 1− m ∏ i=1 e−λit=1−e−Λt. Concluimos pois que a estrutura do relógio adoptada implica que a sobreposição de mprocessos de Poisson independentes de parâmetros λi,i=1,2,...,m, dão origem a um processo de Poisson de parâmetro Λ= m ∑ i=1 λi. Este facto é aplicado de imediato para o caso do Autómato Temporizado Estocástico em que os eventos passíveis de ocorrerem num dado estado xsão dados pelo mapa x→Γ(x). Basta agora considerar Λ(x) = ∑ ei∈Γ(x) λi, e, neste contexto, temos que a função distribuição de probabilidade da variável tempo entre eventos dependente do estado é dada por G(t,x) = P[Y∗(x)≤t] = 1−e−Λ(x)t. Conclusão idêntica relativamente à variável tempo de vida residual se obtém com os mesmos argumentos, ou seja, P[Y∗(x)≤z+t|Y∗(x)>z] = P[Y∗(x)≤t] = G(t,x). 22 Processos de Markov Finalmente, temos todos os ingredientes necessários para caracterizar as probabilidades de transição de estado numa Cadeia de Markov partindo das de um Autómato Temporizado Estocástico em que a ocorrência dos eventos satisfazem as propriedades de Processos de Poisson, ou seja, calcular a probabilidade de transição do estado xpara o estado x0, ou seja, p(x0,x) = P[X0=x0|X= x]. Como P[X0=x0|X=x] = ∑ ei∈Γ(x) P[X0=x0|E0=ei,X=x]P[E0=ei|X=x], e a probabilidade P[X0=x0|E0=ei,X=x] = p(x0|x,ei)é especificada pelo Autómato Temporizado Estocástico, resta-nos calcular P[E0=ei|X=x]. Seja W=min ej∈Γ(x),j6=i{Yj}. Logo, para i=1,...,m, P[E0=ei|X=x] = P[Yi≤W] = Z∞ 0P[Yi≤W|W=y]fW(y)dy, onde fW(y)é a função de densidade de probabilidade associada à distribuição P[Yi≤W|W=y]. Como P[Yi≤W|W=y] = 1−e−Λi(x)y, com Λi(x) = Λ(x)−λi, tem-se que fW(y) = Λi(x)e−Λi(x)y. Logo, P[Yi≤W] = Λi(x)Z∞ 0e−Λi(x)y−e−Λ(x)ydy =λi Λ(x). Observando que o evento [Yi≤W]coincide com [E0=ei]quando o estado é x, conclui-se que p(ei,x) = λi Λ(x).(3.1) Juntando estes elementos, obtemos que p(x0,x) = ∑ ei∈Γ(x) p(x0|x,ei)λi Λ(x). 3.2 Cadeias de Markov Para formular o modelo de uma Cadeia de Markov, teremos de definir um espaço de estados X, a probabilidade do estado inicial p0(x) = P[X0=x]para todo o x∈X, e as probabilidades de transição de estado p(x0,x)em que xé o estado corrente e x0é o estado seguinte. Os desenvolvimentos nesta secção, foram em parte extraídos de [14,1,4]. Para facilitar a notação, considerem-se os estados enumerados e denote-se por pi j(k)a probabilidade da transição de um passo do estado ipara o estado jna k-ésima ocorrência de eventos, ou seja, pi j(k) = P[Xk+1=j|Xk=i]. Obviamente, tem-se que pi j(k)∈[0,1]e que ∑ ∀j pi j(k) = 1. A probabilidade de uma transição de npassos é dada por pi j(k,k+n) = P[Xk+n=j|Xk=i]. Observe-se que, usando a regra da probabilidade total, ou seja, entrando com todas as possibilidades de se passar do estado ipara o estado jpor qualquer estado intermédio rnuma dada 3.2 Cadeias de Markov 23 ocorrência intermédia u(i.e., k<u≤k+n, tem-se que pi j(k,k+n) = ∑ ∀r P[Xk+n=j|Xu=r,Xk=i]P[Xu=r|Xk=i] =∑ ∀r P[Xk+n=j|Xu=r]P[Xu=r|Xk=i], sendo esta última igualdade validada pela propriedade de ausência de memória. Assim, podemos escrever, a conhecida equação de Chapman-Kolmogorov na forma pi j(k,k+n) = ∑ ∀r pir(k,u)pr j(u,k+n). A Matriz de Transição de Estado H(k,k+n) = [pi j(k,k+n)],i,j=0,1,2,..., sendo que irepresenta o índice da linha e jo da coluna, não é mais do que a matriz das probabilidades de transição entre qualquer par de estados, e a equação de Chapman-Kolmogorov toma a forma: H(k,k+n) = H(k,u)H(u,k+n). Uma Cadeia de Markov diz-se homogénea se a probabilidade de transição não depender de k, ou seja, P[Xk+1=j|Xk=i] = pi j. Exemplo. Vamos agora considerar um simples exemplo de um Sistema de Manufactura que iremos modelizar por uma Cadeia de Markov. O sistema de manufactura constituido por duas máquinas - M1 e M2 - com características idênticas para processar tarefas, e operando de acordo com as seguintes regras: •Apenas uma tarefa pode ser submetida no sistema em cada período. •A tarefa pode ser executada por qualquer processador que estiver livre. Se ambas estiverem livres, a tarefa é executada por M1. Se M1 e M2 estiverem ocupadas, a tarefa perde-se. •Se uma tarefa chega e M1 e M2 estão ambos ocupados mas, pelo menos uma delas completa a tarefa que tiver no período de chegada, então a tarefa será processada. Seja pa probabilidade de uma tarefa aparecer no sistema, e qa de uma máquina ocupada terminar a tarefa num qualquer período. Sendo o estado determinado pelo número de máquinas ocupadas, claramente, tem-se que X= {A,B,C}, onde Aé o estado em que M1 e M2 estão livres, Bé o estado em que apenas uma das máquinas está ocupada, e Cé o estado em que ambas as máquinas estão ocupadas. Na figura 3.1, representa-se um autómato em que se descreve a operação do sistema Considerando os diversos eventos que levam à transição para os diversos estados, não é difícel calcular as probabilidades de transição entre os diversos estados, que são dadas por: •pAA =1−p,pAB =p,pAC =0 •pBA = (1−p)q,pBB =pq +(1−p)(1−q),pBC =p(1−q) 24 Processos de Markov Figure 3.1: Exemplo de Sistema de Manufactura modelizado por Cadeia de Markov •pCA = (1−p)q2,pCB =2q(1−q)(1−p)+ pq2,pCC = (1−q)2+2pq(1−q). Neste caso a matriz de transição de estado é dada por P=   pAA pAB pAC pBA pBB pBC pCA pCB pCC   . Uma propriedade importante para Cadeias de Markov diz respeito à distribuição de probabilidade de permanência num dado estado. Sendo V(i)o número de vezes consecutivas que sistema permanece no estado i, ou seja, [V(i) = n] = [Xk+1=i,Xk+2=i,...,Xk+n−1=i,Xk+16=i|Xk=i]. Para se concluir que P[V(i) = n] = (1−pii)pn−1 ii , basta observar que P[V(i) = n] = P[Xk+1=i,Xk+2=i,...,Xk+n−1=i,Xk+n6=i|Xk=i] =P[Xk+n6=i|Xk+n−1=i,...,Xk=i]P[Xk+n−1=i,...,Xk+1=i|Xk=i] = (1−P[Xk+n6=i|Xk+n−1=i,...,Xk=i])P[Xk+n−1=i,...,Xk+1=i|Xk=i] = (1−pii)pn−1 ii , onde o primeiro factor resulta de imediato da propriedade de ausência de memória, e o segundo resulta da conjugação desta propriedade com a aplicação recursiva da regra da probabilidade condicionada: P[Xk+n−1=i,...,Xk+16=i|Xk=i] = P[Xk+n−1=i|Xk+n−2=i,...,Xk=i] ·P[Xk+n−2=i,...,Xk+1=i|Xk=i] =pii ·P[Xk+n−2=i,...,Xk+1=i|Xk=i]. 3.3 Resposta de Cadeias de Markov 25 3.3 Resposta de Cadeias de Markov A resposta de uma Cadeia de Markov consiste na evolução da probabilidade de cada estado à medida que os eventos ocorrem. Seja k=0,1,2,... o contador de ocorrência de eventos e π(k) = [π0(k)π1(k)π2(k)... ]o vector (linha) das probabilidades dos estados na ocorrência k. Nesta secção, consideramos apenas Cadeias de Markov estacionárias, ou seja aquelas em que a Matriz de Transição de Estado não depende da variável de ocorrência k, ou seja, H(k,k+1) = P, para uma matriz Pque caracterização as transições de estado (tal como no exemplo da secção anterior). Dado o vector π(0), a evolução do sistema é dada pela evolução do vector de probabilidades π(k)que, de acordo com as Equação de Chapman-Kolmogorov, satisfaz, para k=0,1,2,..., a equação recursiva π(k+1) = π(k)P.(3.2) Esta equação não é mais do que a versão vectorial da equação que cada uma das suas componentes satisfaz. Cada uma destas é obtida condicionando o evento [Xk+1=j]nos eventos [Xk=i]para todos os valores de i. πj(k+1) = P[Xk+1=j] = ∑ ∀i P[Xk+1=j|Xk=i]P[Xk=i] = ∑ ∀i pi jπi(k). Reiterando a equação (3.2), chega-se à função resposta (em função de k=0,1,2,...) da Cadeia de Markov que é dada por π(k) = π(0)Pk. Esta equação é a base na qual assenta a análise da resposta - quer em regime transitório que em regime permanente - da Cadeia de Markov. Para sistematizar esta análise, importa definir classes de estados apropriadas, o que, por sua vez, requer um conjunto de definições, [1]. Consideraremos apenas as seguintes: •Atingibilidade - o estado jé atingível a partir do estado ise ∃n≥1 tal que pn i j >0, (pn i j é a probabilidade de transitar do estado ipara o estado jem npassos). •Conjunto de Estados Fechado - O conjunto de estados F⊂Xé fechado se pi j =0, ∀(i,j)∈ (F×X\F). •Estado Absorvente - Um estado ié absorvente se {i}é um conjunto fechado. Obviamente, esta proriedade equivale a que pii =1. •Conjunto de Estados Irredutível - Um conjunto de estados fechado Fé irredutível se ∀i,j∈ F,ié atingível de j. Se Xé irredutível, então a Cadeia de Markov é irredutível. •Estados Transitório eRecorrente. Um estado ié dito transitório se ρi<1 e transitório se ρi=1, onde ρié a probabilidade de voltar ao estado idepois de ter lá estado. Sendo ρk i=P[Tii =k], onde a variável tempo de atingimento Ti j é o menor intervalo de tempo para 26 Processos de Markov transitar de ipara j, ou seja, Ti j =min X0=i,Xk=j k>0{k}, tem-se que ρi= ∞ ∑ k=1 ρk i. Observe-se que ρinão é mais do que P[Tii <∞], designando-se Tii de variável tempo de recorrência. •Estados Positivamente Recorrente e de Recorrência Nula. Um estado é de recorrência positiva ou nula conforme Mi<∞ou Mi=∞, onde Mi=E[Tii] = ∞ ∑ k=1 kρk ié o valor esperado do tempo de recorrência. •Estados Periódico eAperiódico. Um estado é periódico se o maior divisor comum do conjunto {n∈IN :pn ii >0}é maior do que 1 e aperiódico doutro modo. É possível demonstrar os seguintes resultados, [1,10]: •Qualquer Cadeia de Markov com um espaço dos estados finito tem pelo menos um estado recorrente. •Se o estado ié recorrente e jé atingível de i, então jé recorrente. Consequentemente, qualquer estado de um conjunto fechado e irredutível é recorrente. •Se o estado ié positivamente recorrente e jé atingível de i, então jé positivamente recorrente. Consequentemente, qualquer estado de um conjunto fechado e irredutível é positivamente recorrente. •Para qualquer conjunto fechado e irredutível, qualquer estado é positivamente recorrente, de recorrência nula, ou transitório. •Se uma Cadeia de Markov é irredutível, então todos oss estados têm o mesmo período. Por resposta em regime permanente, entende-se na determinação nas probabilidades dos diversos estados quando k→∞, ou seja, para o estado i,πj=lim k→∞πj(k). Caso este limite exista e verifique ∑ ∀i πi=1, então o vector π= (π0,π2,...)designa-se de regime permanente ou ponto de equilíbrio. Extraindo o limite k→∞na equação de transição de estado π(k+1) = π(k)P, obtém-se π=πP, onde as suas componentes têm que se não negativas e ∑ ∀i πi=1. A análise do regime permanente é considerada em duas situações: a) Cadeias de Markov irredutíveis. Importa observar que se existirem estados periódicos numa Cadeia de Markov irredutível, então, obviamente o limite de π(k)quando k→∞não está definido. É possível demonstrar o seguinte resultado: Teorema 1. Para qualquer Cadeia de Markov irredutível e aperiódica, ∀j, o limite lim k→∞πj(k)existe e é igual a um valor πjque não depende de πj(0). 3.3 Resposta de Cadeias de Markov 27 Se todos os estados desta Cadeia de Markov forem transitórios ou de recorrência nula, então πj=0 e, portanto, o limite πnão constitui uma probabilidade estacionária. Se todos os estados desta Cadeia de Markov forem positivamente recorrentes (também referidos pro estados ergódicos), então πj=1 Mj . Neste caso, o limite πexiste e constitui uma probabilidade estacionária. Considere agora a seguinte cadeia de nascimento-morte representada no diagrama 3.2. Figure 3.2: Cadeia Nascimento-Morte Seja a probabilidade ptal que 0 <p<1, então resolvendo a equação π=πP, tem-se que: π0=π0p+π1pou seja π1=1−p pπ0 πj=πj−1(1−p)+πj+1p,j=1,2,.... Repare que se tem π1=π0(1−p) + π2po que determina que π2=1−p p2 π0. Resolvendo por indução, chega-se à conclusão que πj=1−p pj π0,j=1,2,.... Utilizando o facto de que ∞ ∑ j=0 πj=π0 ∞ ∑ j=01−p pj =1, conclui-se que π0= ∞ ∑ j=01−p pj!−1 . Vamos agora considerar as diversas possibilidades para os valores de p: i) 1−p p<1 ou seja p>0.5. Neste caso, ∞ ∑ j=01−p pj =p 2p−1e πj=2p−1 p1−p pj ,j=0,1,2,.... ii) Se p<0.5, tem-se que a série ∞ ∑ j=01−p pj diverge e, portanto, πj=0, j=0,1,2,.... Neste caso a cadeia de Markov é transitória. iii) Se p→0.5, conclui-se que o estado 0 é de recorrência nula. b) Cadeias de Markov redutíveis. Neste caso, duas possibilidades podem ser consideradas: existe apenas um subconjunto de estados irredutíveis, ou, então existem vários destes conjuntos. No primeiro caso, o estado entra, mais tarde ou mais cedo, no conjunto irredutível e 28 Processos de Markov fica lá daí em diante, reduzindo-se à situação analisada em a). No segundo, coloca-se apenas a questão de saber qual é o subconjunto irredutível no qual o estado vai entrar primeiro, pois o estado nele permanecerá daí em diante. De qualquer modo, verifica-se que basta apenas determinar qual a probabilidade do estado entrar num dado subconjunto de estados irredutíveis (que designamos de I), a partir de um dado estado inicial no subconjunto de estados transitórios (que designamos de T). Seja ρi(I) = P[Xi∈I|X0=i]. Utilizando probabilidades condicionadas e a propriedade de ausência de memória, tem-se que: ρi(I) = ∑ j∈I Pi j +∑ l∈T ρl(I)Pil. Em determinadas condições, nomeadamente, para subconjuntos Tfinitos, é possível resolver este sistema de equações e determinar ρi(I). 3.4 Cadeias de Markov em Tempo Contínuo Nesta secção vamos considerar a situação em que as probabilidades de transição da Cadeia de Markov variam em tempo contínuo em vez de por passos discretos em que o tempo é irrelevante, ou seja, agora tem-se que, para s≤t,pi j(s,t) = P[X(t) = j|X(s) = i]. Neste contexto, a propriedade de Markov formulada para o contexto discreto (2.3) toma agora a forma: P[X(tk+1) = xk+1|X(tk) = xk,X(tk−1) = xk−1,...,X(t0) = x0] =P[X(tk+1) = xk+1|X(tk) = xk].(3.3) Embora a análise tenha analogias com o caso discreto, agora não é mais possível considerar a matriz de transição de um passo P, uma vez que as probabilidades de transição de estado vão depender do tempo que tiver decorrido. Para deduzir a equação de Chapman-Kolmogorov neste caso, adapta-se a metodologia no caso anterior, obtendo-se, pela regra da probabilidade total e pela propriedade de ausência de memória, pi j(s,t) = ∑ ∀r P[X(t) = j|X(u) = r,X(s) = i]P[X(u) = r|X(s) = i] =P[X(t) = j|X(u) = r]P[X(u) = r|X(s) = i] =∑ ∀r pir(s,u)pr j(u,t) Logo, designando a matriz [pi j(s,t):∀i,j]por H(s,t), tem-se, para s≤u≤t, que H(s,t) = H(s,u)H(u,t). Importa determinar como é que a matriz de transição de estado varia com o tempo, ou seja, calcular 3.4 Cadeias de Markov em Tempo Contínuo 29 ∂H(s,t) ∂t. Note-se que, para s≤t, ∂H(s,t) ∂t=lim δ→0 H(s,t+δ)−H(s,t) δ =lim δ→0 H(s,t)H(t,t+δ)−H(s,t) δ =H(s,t)lim δ→0 H(t,t+δ)−I δ =H(s,t)Q(t) onde Q(t) = lim δ→0 H(t,t+δ)−I δé a matriz taxa de transição de probabilidades. Resolvendo esta última equação diferencial, obtém-se H(s,t) = eRt sQ(τ)dτ. Tal como no caso discreto, vamos apenas considerar apenas Cadeias de Markov homogéneas. Neste caso, tem-se que a probabilidade de transição pi j(s,t)depende apenas de t−s, ou seja, pi j(τ) = P[X(s+τ) = j|X(s) = i]e, analogamente ao caso discreto, considera-se P(τ) = H(s,s+ τ). Observe-se que, agora a matriz taxa de transição de probabilidades não depende de te é, portanto, constante, ou seja, Q. Logo, para o caso homogéneo, a equação de Chapman-Kolmogorov na forma diferencial em tempo contínuo é dP(τ) dτ=P(τ)Q,(3.4) sendo a sua solução dada por P(τ) = eQτ.(3.5) Utilizando argumentos similares aos do caso discreto, conclui-se que a função distribuição de probabilidade da variável aleatória tempo de permanência num estado é exponencial, ou seja, para t≥0, é dada por P[V(i)≤t] = 1−e−Λ(i)t. onde Λ(i) = ∑ ei j∈Γ(i) λi j, sendo λi j a taxa de Poisson associada ao evento ei j viável no estado i. Tal como no caso discreto, a propriedade de ausência de memória - P[V(i)>s+t|V(i)>s] = P[V(i)>t]- também se verifica. Importa interpretar cuidadosamente a matriz Qpara clarificar as equações (3.4) e, por consequência, a sua solução dada por (3.5). Sendo Q= [qi j], a equação (3.4) pode ser escrita por extenso na forma dpii(τ) dτ=pii(τ)qii +∑ r6=i pir(τ)qri. Fazendo τ=0, tem-se que −qii =−d pii(τ) dτ|τ=0=d[1−pii(τ)] dτ|τ=0. Uma vez que 1 −pii(τ) 36 Teoria das Filas de Espera a) Utilização do sistema - Fracção do tempo em que os servidores estão ocupados. No caso de filas com apenas um servidor, tem-se que, em regime permanente, Utilização do sistema =1−π0. b) Fluxo do sistema - Taxa a que os clientes saiem do sistema. Observe-se que, no caso de filas com apenas um servidor, tem-se que, em regime permanente, Fluxo do sistema =µ(1−π0). c) Intensidade de tráfego - Fracção da taxa média de chegada sobre a taxa média de serviço. Esta medida de desempenho designa-se por ρe, para o caso Markoviano considerado anteriormente, têm-se ρ=λ µ. No caso de existirem qservidores, o denominador será dado por µ1+···+µq. Observe-se que, no caso de filas com apenas um servidor, o sistema só está em regime permanente se λ=µ(1−π0), donde se conclui que, para estes sistemas a Intensidade de Tráfego ρé igual a 1−π0. As variáveis acima indicadas constituem a base para a definição das equações genéricas da dinâmica da fila de espera. Uma vez que Dk−1≤Akimplica que Wk=0, tem-se que Wk= max{0,Dk−1−Ak}. Por outro lado, como Yk=Ak−Ak−1, pode-se escrever a equação da dinâmica relativa à variável tempo de espera que é dada por Wk=max{0,Wk−1+Zk−1−Ak}.(4.1) De forma semelhante, usando as definições anteriores, tem-se que as variáveis tempo no sistema e instante de partida são dadas por: Sk=max{0,Sk−1−Yk}+Zk(4.2) Dk=max{Ak,Dk−1}+Zk(4.3) Um resultado importante para filas de espera arbitrárias em regime permanente é a chamada Lei de Little, [18], E[X] = λE[S]. Para o plausibilificar, considere-se o espaço de eventos E={a,d}, sendo a variável aleatória (VA) comprimento da fila X(t), dada por X(t) = Na(t)−Nd(t), onde NaeNdsão os processos de contagem associados aos eventos aedno intervalo (0,t], respectivamente. Consideremos agora uma realização específica destas variáveis (relembra-se que as letras maiúsculas são reservadas para as VA e as minúsculas para os seus valores), onde u(t)é o valor do tempo total dos clientes no sistema. Facilmente se conclui que, para este intervalo, o tempo médio de cada cliente no sistema, o comprimento médio da fila, e a taxa média de chegada de clientes, são, para esta realização, dados, respectivamente, por ¯s(t) = u(t) na(t), ¯x(t) = u(t) t, e λ(t) = na(t) t. 4.3 Fila de espera M/M/1 37 Daqui se conclui que ¯x(t) = λ(t)¯s(t). Assumindo agora que, para esta realização, se tem que lim t→∞λ(t) = λe que lim t→∞¯s(t) = ¯s, conclui-se que ¯x=λ¯s. Se o sistema for ergódico, ou seja, se estes limites existirem para qualquer realização, conclui-se que E[X] = λE[S]. Argumentos de idêntica natureza também plausibilificam E[XQ] = λE[W]eE[XS] = λE[Z], onde XQeXS, são as VAs descrevendo os elementos da fila que, respectivamente, esperam para ser servidos e estão a ser servidos. 4.3 Fila de espera M/M/1 Nesta secção, apresenta-se a dedução das medidas de desempenho de uma fila simples do tipo M/M/1, ou seja, os processos de entrada e de serviço são Markovianos, só existe um servidor e a capacidade da fila é ilimitada. O espaço dos eventos e dos estados são, respectivamente E={a,d}eX={0,1,2,...}, representando a variável aleatória X(t)o comprimento da fila no instante te{X(t)}é um processo estocástico Generalizado de Semi-Markov. Sendo muito difícil obter a caracterização do regime estacionário para sistemas tão gerais, iremos apenas considerar filas em que os processos de chegada e de partida são Markovianos, ou seja, as funções distribuição de probabilidade dos tempos entre chegadas e de entre serviços dados, respectivamente por, A(t) = 1−e−λte B(t) = 1−e−µt, onde λeµsão parâmetros dados. Importa observar que, para x>0, se tem (f(x,a) = x+1 f(x,d) = x−1 podemos considerar esta fila de espera como uma cadeia nascimento-morte, tendo já sido concluido em (3.9) que πn=λ0...λn−1 µ1...µnπ0,n=1,2,... (4.4) π0=1 1+∑∞ n=1λ0...λn−1 µ1...µn(4.5) Nesta notação, não aparece qualquer indicação relativa à disciplina de fila e, neste caso, geralmente, assume-se a política de “O primeiro a chegar é o primeiro a ser servido”. No caso de uma fila do tipo M/M/1, existe apenas um servidor, a capacidade da fila é ilimitada, e os processos de chegada e de partida são Markovianos, sendo as funções distribuição de probabilidade dos tempos entre chegadas e de entre serviços dados, respectivamente por, A(t) = 1−e−λt eB(t) = 1−e−µt, onde λeµsão parâmetros dados. No que se segue, vai ser importante o facto de que, para processos de chegada de Poisson, se tem que P[chega cliente em te o estado é X(t) = n] = P[o estado do sistema em téX(t) = n]. 38 Teoria das Filas de Espera Vamos agora deduzir as fórmulas da resposta do sistema e das mais importantes funções de desempenho para a fila do tipo M/M/1 com λj=λeµj=µ,∀j, cujo diagrama está representado na figura 4.2. Figure 4.2: Transição de estados para fila M/M/1 Neste caso, e exprimindo em termos do indicador Intensidade de Tráfego ρ=λ µdefinido anteriormente nesta secção, (4.4)e(4.5) tomam, respectivamente, as formas πn=ρnπ0,n=1,2,... (4.6) π0=1 1+∑∞ n=1ρn(4.7) Caso ρ<1, ou seja λ<µ, a série em (4.7) converge, obtendo-se ∞ ∑ n=1 ρn=ρ 1−ρ. Donde se conclui que πn=ρn(1−ρ),n=0,1,2,.... As expressões para as funções de desempenho em regime estacionário tipicamente consideradas e referidas na secção anterior, são obtidas directamente, ou extraídas de forma imediata utilizando resultados convencionais de convergência de séries. Assim: •Intensidade de tráfego -ρ=λ µ •Utilização do sistema - 1−π0=ρ •Fuxo do sistema -lambda =µ(1−π0) •Comprimento da fila médio -E[X] = ∞ ∑ n=1 nπn=ρ 1−ρ. Para concluir esta expressão, observe-se que ∞ ∑ n=1 nπn= (1−ρ) ∞ ∑ n=1 nρn= (1−ρ)ρ ∞ ∑ n=0 dρn dρ = (1−ρ)ρd dρ ∞ ∑ n=0 ρn!= (1−ρ)ρd(1−ρ)−1 dρ = (1−ρ)ρ(1−ρ)−2. 4.3 Fila de espera M/M/1 39 •Tempo no sistema médio -E[S] = 1 λE[X] = 1 µ(1−ρ). •Tempo de espera médio -E[W] = ρ µ(1−ρ). Para concluir esta expressão observe-se que, sendo E[XQ] = E[X]−E[XS],E[W] = 1 λE[XQ], e que E[XS] = λE[Z], tem-se que E[W] = 1 λE[X]−E[Z] = ρ λ(1−ρ)−1 µ=ρ µ(1−ρ). Argumentos semelhantes, permitem obter as funções de resposta, ou seja, as probabilidades dos estados, e as funções de desempenhos para os outros tipos de fila mencionados neste capítulo. Contudo, para efeitos desta dissertação, vamos restringir-nos à filas da classe M/M/1. 40 Teoria das Filas de Espera Chapter 5 Aplicação ao Controlo e Optimização de Sistemas de Manufactura Discretos Neste capítulo, apresenta-se a abordagem para aplicar os conceitos e resultados anteriores para modelizar e optimizar o regime permanente de Sistemas de Manufactura. Apesar dos conceitos e resultados apresentados nos capítulos anteriores serem muito gerais, vamos, como já foi referido, apenas considerar filas de espera do tipo M/M/1. Assim, numa primeira secção, apresenta-se a modelização de um sistema de produção discreto como uma rede de filas de espera. No sentido de limitar a complexidade para não obscurecer o objectivo essencial desta disssertação, vamos apenas considerar redes de filas de espera abertas, ou seja, aquelas para as quais o número de clientes no sistema não tem que permanecer necessáriamente constante. Esta secção vai apresentar um exemplo que será tratado no decorrer do capítulo. Na secção seguinte, apresentar-se-ão alguns resultados específicos para redes de filas de espera que constituem uma extensão resultados abordados nos capítulos anteriores. Finalmente, na última secção, vamos tratar o exemplo introduzido previamente e obter uma solução que optimiza um critério de desempenho seleccionado. 5.1 Modelização e Análise de um Sistema de Manufactura Um sistema de produção discreto consiste num conjunto de máquinas discretas interligadas nas quais ocorrem operações sobre componentes do produto ou produtos a fabricar. Tipicamente, cada máquina pode ser vista como uma fila de espera constituída por um ou mais servidores que executam uma operação, bem como por um inventário à entrada (referido como “fila” no senso comum) onde as peças esperam para ser trabalhadas. Uma vez que o custo de produção aumenta fortemente com o tempo de permanência das peças no sistema, importa que este inventário distribuído seja o menor possível. Contudo, o carácter aleatório dos processos de produção que se manifesta de várias formas - duração das operações, taxa de sucesso destas que tem de ser verificada pelo controlo de qualidade, avarias, etc - torna importante a existência do inventário distribuído por forma 41 42 Aplicação ao Controlo e Optimização de Sistemas de Manufactura Discretos a minimizar a ruptura de stocks e assim contribuir para uma maior robustez do funcionamento do sistema. Assim, a especificação de um sistema de manufactura requer os seguintes ingredientes: •Taxas de entrada de clientes do exterior no sistema. Estes podem ser peças ou unidades de matéria prima a serem trabalhadas nas máquinas que os vão processar. •Cada máquina é representada por uma fila de espera, ou seja, por um servidor que realiza as operações de produção a uma dada taxa e por um “armazém” à entrada no qual as peças aguardam a disponibilidade do servidor. •Taxa de serviço do servidor de cada uma das máquinas. No caso destas estarem modelizadas por processos de Poisson, esta taxa não é mais do que o parâmetro que caracteriza a respectiva distribuição probabilística. •Probabilidade de transição do cliente uma vez processado numa dada máquina para a máquina seguinte, ou, eventualmente, para a saída do sistema. As transições entre máquinas são especificadas pelo “routing” que define o sistema de produção. Observe que, no caso da saída de uma máquina tiver como destino apenas uma outra máquina, esta probabilidade será necessáriamente 1. Por outro lado, a soma das probabilidades de todas as transições possíveis da saída de um dada máquina terá quer ser igual a 1. Importa referir que uma dada operação de produção pode não preencher os requisitos pretendidos e necessitar de ser repetida com uma certa probabilidade. •Especificação da saída ou saídas do sistema, através das quais os clientes deixam o sistema, normalmente correspondendo a peças acabadas ou, então, a peças rejeitadas pelo sistema no sentido de não ser possível qualquer processamento adicional. •Variáveis de “controlo” que poderão ser as taxas das diversas entradas ou então os parâmetros que definem a taxa de transição entre máquinas. O valor destes parâmetros poderão ser escolhidos por forma a optimizar uma ou mais funções de desempenho do sistema. A título de exemplo considere-se uma sistema de produção discreto representado no diagrama 5.1. Observe que, neste sistema existem 4 máquinas e que as probabilidades de transição das peças entre estas têm que satisfazer as seguintes relações: pi j ∈[0,1],p14 +p13 =1, p21 +p22 +p23 =1 ep33 +p34 =1. Repare que existem duas entradas que poderão corresponder a peças ou matérias primas diferentes ou que requerem operações iniciais diferentes e que existe apenas uma saída após uma operação final na máquina M4. Finalmente, importa observar que as máquinas M2, M3 e M4 estão realimentadas. Utilizando os resultados abordados nos capítulos anteriores desta dissertação, o tratamento deste exemplo envolve as seguintes componentes: 5.2 Redes de Filas de Espera Markovianas 43 Figure 5.1: Sistema de produção discreto representado por uma rede de filas de espera. •Determinação da taxa de chegada total de peças a cada máquina, λi,i=1,2,3,4, por forma a que o sistema esteja em equilíbrio em regime estacionário. Estas taxas são determinadas em função dos parâmetros do sistema acima indicados. •Escolha do valor dos parâmetros livres por forma a optimizar uma dada função de desempenho. Este exemplo irá ser tratado na última secção deste capítulo. Contudo, para o efeito, precisamos ainda de alguns resultados que estendem os dos capítulos anteriores para tratar redes de fila de espera. 5.2 Redes de Filas de Espera Markovianas Nesta secção, começamos por caracterizar o processo de partida de uma fila do tipo M/M/1, [12]. Considere-se a variável intervalo de tempo entre as partidas kek−1, Ψk=Dk−Dk−1, para K=1,2,..., com Ψ0=0, e assume-se que existe a correspondente distribuição de probabilidade em regime estacionário P[Ψ≤t] = lim k→∞P[Ψk≤t]. Neste contexto, vamos determinar P[Ψ≤t]. De (4.3), conclui-se que Ψk=Dk−Dk−1=max{Ak−Dk−1,0}+Zk. Logo, Ψk=(Zkse Ak≤Dk−1 Ak−Dk−1+Zkse Ak>Dk−1 Repare-se que, no primeiro caso, o instante de partida do cliente depende apenas do servidor que está sempre ocupado e que esta situação se verifica quando a “fila" não está vazia no instante imediatamente anterior à k-ésima chegada, ou seja X(A− k)>0. No segundo caso, o servidor está desocupado, o que se verifica quando X(A− k) = 0, e o instante de partida é dado pela soma do tempo de serviço e a diferença entre os instantes de chegada do cliente corrente e de saída do 44 Aplicação ao Controlo e Optimização de Sistemas de Manufactura Discretos cliente anterior. Assim, tem-se que P[Ψ≤t] = P[Ψk≤t|X(A− k)>0]P[X(A− k)>0]+P[Ψk≤t|X(A− k) = 0]P[X(A− k) = 0] =P[Zk≤t]P[X(A− k)>0]+P[Ak−Dk−1+Zk≤t]P[X(A− k) = 0]. Uma vez que a fila é Markoviana temos que P[Zk≤t] = 1−e−µt. Por lado, como estamos a considerar apenas o regime permanente, tem-se lim k→∞P[X(A− k) = 0] = π0=1−ρe, portanto, lim k→∞P[X(A− k)>0] = ρ. Finalmente, resta calcular a função distribuição de probabilidade da soma das VAs Ak−Dk−1eZk. Sendo a função densidade de probabilidade correspondente dada pela convolução das funções densidade de probabilidade das VAs parcelas e, como se tem P[Zk≤t] = 1−e−µte, pela propriedade de Markov, P[Ak−Dk−1≤t] = 1−e−λt, obtém-se P[Ak−Dk−1+Zk≤ t] = λe−µt−µe−λt µ−λ. Combinando as probabilidades, conclui-se que, para uma fila M/M/1 em regime estacionário, Ψtem uma distribuição de Poisson com parâmetro λ, ou seja, P[Ψ≤t] = 1−e−λt. Este resultado é fundamental para se determinar a taxa total de chegada para cada um dos nós de uma rede em regime estacionário e assim escrever as equações desta. Consideremos agora uma rede com Nfilas de espera do tipo M/M/1 e subordinadas à disciplina “primeiro chegado, primeiro servido”, sendo cada uma vista como um nó. Sejam •λieµi,i=1,...,N, os parâmetros de Poisson dos respectivos processos de chegada e de serviço, •pi j, a probabibilidade de um cliente acabado de ser servido num dado nó itransitar para o nó j, e •ri, a probabilidade de chegada de um cliente do exterior do sistema ao nó i. A taxa total de chegada no nó i,i=1,...,N, é dada por λi=ri+ N ∑ j=1 λjpji. Observe que, nesta formulação, existe a possibilidade de “feedback”, ou seja, de um cliente acabado de ser servido no nó ivoltar a entrar neste com uma certa probabilidade. Importa também referir que se demonstra que a estabilidade é garantida se λi<µi,i=1,...,N. Importa observar que, embora a sobreposição do processo de chegada externo com o processo de saída realimentado não seja de Poisson, demonstra-se que o processo de partida é de Poisson e com taxa igual à do processo de chegada externo. Consideremos agora o caso de uma rede formada por duas filas de espera em série representadas na figura 5.2. 5.2 Redes de Filas de Espera Markovianas 45 Figure 5.2: Diagrama para duas filas M/M/1 em série Seguindo o formalismo de [1], consideramos dois nós cuja variável comprimento de fila, [X1,X2], tem dimensão 2, sendo o espaço dos estados X={(n1,n2):n1=0,1,2,...,n2=0,1,2,...} e o espaço de eventos E={a,d1,d2}. Considerando λ,µ1eµ2as taxas dos processos de Poisson, respectivamente, de chegada e de serviço de cada um dos dois servidores, podemos escrever o sistema de equações:                0=λπ(n1−1,n2)+ µ1π(n1+1,n2−1)+ µ2π(n1,n2+1) −(λ+µ1+µ2)π(n1−1,n2)se n1>0,n2>0 0=λπ(n1−1,0)+ µ2π(n1,1)−(λ+µ1)π(n1−1,0)se n1>0,n2=0 0=µ1π(1,n2−1)+ µ2π(0,n2+1)−(λ+µ2)π(0,n2)se n1=0,n2>0 0=µ2π(0,1)−λπ(0,0)se n1=n2=0. Dada a independência dos processos, temos que π(n1,n2) = π(n1)π(n2)e o resultado anterior respeitante à caracterização do processo de partida, as probabilidades na equação podem ser escritas em termos das dos processos nascimento-morte de uma dimensão, ou seja, sendo ρi=λ µi para i=1,2, tem-se que π(n1,n2) = (1−ρ1)ρn1 1(1−ρ2)ρn2 2. A generalização para Nnós é imediata, bastando observar que agora se tem π(n1,n2,...,nN) = π1(n1)π2(n2)···πN(nN). Temos já todos os ingredientes para determinar a probabilidade π(n1,n2,...,nN)em regime estacionário de para dada rede de filas de espera do tipo M/M/1 com a disciplina “primeiro chegado, primeiro servido”. O método envolve os seguintes passos: 1) Escrever as equações dos nós. Estas são dadas pelo sistema de Nequações λi=ri+ N ∑ j=1 λjpji, para j=1,2,...,N. 2) Resolver o sistema de equações em ordem aos λi,i=1,2,...,N, que serão as taxas de chegada em cada nó com o sistema em equilíbrio em regime estacionário. 3) Para cada i, calcular o parâmetro ρi. 4) Calcular π(n1,n2,...,nN). 52 REFERENCES [17] L. Kleinrock. Queueing Systems. Wiley, New York, 1975. [18] J. Little. A proof of L=λW.Operations Research, 9 (3):383–387, 1961.