Optimizacija mravlje kolonije

Sad prelazimo sa analogije metalurških procesa na biologiju :) E ta veštačka inteligencija, uvek traže inspiraciju u nečemu postojećem. Ideja je nastala na osnovu proučavanja konunikacije i upravljanja koji se dešava među insektima. Mravi kad traže neki izvor ispuštaju feromone. Feromoni vremenom isparavaju a količina feromona na nekoj stazi zavisi i od dužine staze. Mrav koji nađe kraći put na tom putu će ostaviti veću količinu feromona koji će drugim mravima poslužiti za navođenje na bolju putanju.

Pomoću ove ideje moguće je napisati heuristiku za rešavanje problema trgovačkog putnika (nalaženje što kraće Hamiltonove konture). 

Dovoljno je imati 3 proste formule, toliko proste da su matematičari morali da ih zakomplikuju lupanjem nekih čudnih oznaka, samo im rusko slovo Я fali. To me podseti na epizodu kad sam polagao statistiku kod jedne nemaštovite asistentkinje kojoj je bilo čudno što formulu uslovne verovatnoće nisam izveo identično kao u knjizi, tj. što nisam bubao već pravio delikt mišljenja. Pa sam joj rekao da dok god je formula ispravna mogu da koristim i Э kao oznaku :D Matematičar su uvek najveći problem matematike.

Zato ću ja formule prikazati uprošteno (a i ne znam kako bih napisao sve one indekse, stepene, klince, palce). Prva formula služi da se izračuna količina feromona za neku ivicu (direktnu vezu dva grada). Formula je rekurentna (programerski rečeno rekurzivna) jer se nova količina feromona izračunava preko stare prostim dodavanjem:

kolicina_feromona = kolicina_feromona + Q/duzina_Hamiltonove_konture *RO

Ili uprošćeno:

kolicina_feromona +=  Q/duzina_Hamiltonove_konture *RO

Q je neka konstanta koja označava količinu feromona kojom svaki mrav raspolaže. Pa mrav koji nađe kraću putanju ima ivicama da doda više feromona jer je Q/mala_vrednost veći broj nego Q/velika_vrednost. Konstanta RO je između 0 i 1 i služi da odredi koliko novog feromona treba dodati. Pošto je formula rekurentna mora postojati početna vrednost pa se uzima da je kolicina_feromona = 1 / broj_gradova

Druga formula simulira isparavanje feromona:

kolicina_feromona = kolicina_feromona * (1 - RO)

ili:

kolicina_feromona *= (1 - RO)

RO je ista konstanta kao iz prethodne formule. Npr. ako je RO = 0,7 to znači da će se na putanju kojom je mrav prošao staviti 70% vrednosti Q/duzina_Hamiltonove_konture, a od stare količine feromona 70% će ispariti.

I treća formula izračunava verovatnoću da se krene nekom ivicom tj. da se obiđe neki neposećeni grad iz tekućeg grada. Za svaku moguću ivicu (tj. za svaku direktnu vezu tekućeg i neposećenih gradova) se izračunava poželjnost. Poželjnost je proporcionalna količini feromona na ivici i njenoj kratkoći. Kratkoću definišemo kao 1 / duzina_ivice a količinu feromona čitamo direktno. Udeo pozeljnosti i-te moguće ivice u zbiru pozeljnosti svih mogućih ivica je upravo verovatnoća da se mrav zaputi tom ivicom. Količinu feromona matematičari označavaju ubercool slovom Tau, a kratkoću slovom Eta. Formula je (** je operator stepenovanja):

P = (Tau**ALFA * Eta**BETA) / SUM(Tau**ALFA * Eta**BETA)

Da bi se kontrolisao značaj kratkoće i feromona uvode se ALFA i BETA parametri kojim se stepenuju količine feromona (Tau) i kratkoće (Eta). Pošto su i Tau i Eta manji od 0 znači da što je stepen veći to se umanjuje značaj te karakteristike u poželjnosti. 

I konačno kod:

import random, copy, math, time

class Grad(object):     
    def __init__(self, x, y, ime):
        self.x = x
        self.y = y
        self.ime = ime
        
    def rastojanje(self, grad):
        """ 
        rastojanje ovog grada od grada primljenog kao prarametar
        preko euklidske metrike (odnosno funkcije hypot).
        """
        return math.hypot(self.x - grad.x, self.y - grad.y)
        
    def __str__(self):
        return "%s(%d, %d)" % (self.ime, self.x, self.y)
 
class Ivica(object):
    def __init__(self, feromoni, duzina):
        self.feromoni = feromoni
        self.duzina = duzina

class Mrav(object):
    # atribut klase u kojem se cuva broj svih stvorenih mrava
    brojac = -1
    gradovi = None
    ivice = None
    
    def inicijalizuj_atribute(self):
        self.tek_grad = self.pocetni_grad
        self.putanja = [self.tek_grad]  
        self.duzina_putanje = 0.0
        
    def __init__(self, gradovi, alfa, beta, ro, q):
        Mrav.ALFA = alfa
        Mrav.BETA = beta
        Mrav.RO = ro
        Mrav.Q = q
        
        if not Mrav.gradovi: 
            Mrav.gradovi = gradovi
            
        if not Mrav.ivice:
            Mrav.ivice = {}
            init_feromoni = 1.0 / len(Mrav.gradovi)
            for i in Mrav.gradovi:
                for j in Mrav.gradovi:
                    if i != j :
                        rastojanje = i.rastojanje(j)
                        Mrav.ivice[i, j] = Ivica(init_feromoni, rastojanje)
                        
        Mrav.brojac += 1
        self.pocetni_grad = Mrav.gradovi[Mrav.brojac % len(Mrav.gradovi)]
        
        self.inicijalizuj_atribute()
       
    def _tau_x_eta(self, do_grad):
        # izracunava vrednost pozeljnosti ivice od tekuceg grada do grada
        # primljenog kao parametar. Pozeljnost ivice je direktno proporcionalna 
        # kolicini feromona na njoj (sto se oznacava sa tau), a obratni proporcionalna 
        # njenoj duzini (eta). Stepenovanjem tau i eta sa parametrima ALFA i BETA 
        # omogucava nam da dajemo prednost jednom ili drugom kriterijumu pozeljnosti
        # i to tako sto BETA oznacava koliko umanjujemo znacaj blizine a ALFA
        # koliko umanjujemo znacaj feromona. Npr. za ALFA: 1 i BETA: 5 znaci da
        # nam je vaznija kolicina feromona od blizine
        ivica = Mrav.ivice[self.tek_grad, do_grad]
        tau = ivica.feromoni
        eta = 1.0 / ivica.duzina
        
        return tau**Mrav.ALFA * eta**Mrav.BETA
        
    def _odaberi_sledeci_grad(self):
        imenilac = 0.0
        
        # svi gradovi razlika poseceni = neposeceni, skupovna operacija
        neposeceni_gradovi = set(Mrav.gradovi).difference(self.putanja)
        
        # ako nema neposecenih vrati None
        if not neposeceni_gradovi:        
            return None
        
        # za svaki neposeceni grad izracunaj _tau_x_eta i saberi
        # tj. racunamo zbir svih pozeljnosti posete neposecenim gradovima
        for grad in neposeceni_gradovi:
            imenilac += self._tau_x_eta(grad)
        
        # za svaki neposecen grad izracunaj verovatnocu posete p
        # i pokusaj da je ostvaris, ako je ne ostvaris probaj
        # ponovo sve iz pocetka dok konacno ne uspes u tome.
        # Verovatnoca posete se izravunava tako sto podelimo pozeljnost posete
        # neposecenog grada sa zbirom svih pozeljnosti posete neposecenim gradovima
        while True:
            for grad in neposeceni_gradovi:
                p = self._tau_x_eta(grad) / imenilac
                
                if random.random() < p:
                    return grad

    def pusti_mrava(self):
        sled_grad = self._odaberi_sledeci_grad()
        
        # dok postoji sledeci grad za obilazak stavi ga u putanju obishenih
        # azuriraj duzinu putanje mrava i tekuci grad a zatim nadji sledeci grad
        while sled_grad:
            self.putanja.append(sled_grad)
            self.duzina_putanje += Mrav.ivice[self.tek_grad, sled_grad].duzina
            self.tek_grad = sled_grad
            sled_grad = self._odaberi_sledeci_grad()
        
        # sad je sled_grad None znaci da smo prosli sve gradove
        # a Hamiltonovu konturu zatvaramo vezujuci 1. i poslednji 
        # poseceni grad koji se cuvaju u self.putanja (0. i -1. poziciju)
        self.duzina_putanje += Mrav.ivice[self.putanja[0], self.putanja[-1]].duzina
        
    def pospi_feromone(self):
        br_gradova = len(self.putanja)
        
        # svakoj ivici putanje mrava (Hamiltonove konture) u oba smera dodaj 
        # odgovarajucu kolicinu feromona po formuli Q / duzina_putanje * RO
        for i in xrange(br_gradova):
            od_grada = self.putanja[i]
            do_grada = self.putanja[(i + 1) % br_gradova]
            
            Mrav.ivice[od_grada, do_grada].feromoni += Mrav.Q / self.duzina_putanje * Mrav.RO
            Mrav.ivice[do_grada, od_grada].feromoni = Mrav.ivice[od_grada, do_grada].feromoni
            
    @staticmethod
    def ispari_feromone():
        # smanji kolicinu feromona na svakoj ivici po rekurentnoj formuli feromoni *= 1 / RO
        for ivica in Mrav.ivice.itervalues():
            ivica.feromoni *= (1.0 - Mrav.RO)
            
    def koordinate(self):
        x, y = [], []
        
        for grad in self.putanja:
            x.append(grad.x)
            y.append(grad.y)
        
        return x, y
    
    # toString u javi        
    def __str__(self):
        return "duzina puta: %f" % self.duzina_putanje
    
    # compare u javi
    def __cmp__(self, drugi):
        return cmp(self.duzina_putanje, drugi.duzina_putanje)
            
            
class MravljiSimulator(object):
    def __init__(self, gradovi, ALFA = 1.0, BETA = 5.0, RO = 0.5, Q = 100.0):
        self.mravi = []
        for i in xrange(len(gradovi)):
            self.mravi.append(Mrav(gradovi, ALFA, BETA, RO, Q))
        self.najbolji_mrav = self.mravi[0] # proizvoljnog mrava proglasimo najboljim
        
    def start(self, ponavljanja=20):
        # za svako ponavljanje
        for i in xrange(ponavljanja):
            # pusti svakog mrava da napravi Hamiltonovu konturu
            for mrav in self.mravi:
                mrav.pusti_mrava()
                
            Mrav.ispari_feromone()
                
            # svaki mrav nek pospe feromone na svojoj konturi
            for mrav in self.mravi:
                mrav.pospi_feromone()
                
            najuspesniji_mrav_generacije = min(self.mravi)
            
            # pamtimo sveukupno najuspesnijeg mrava
            self.najbolji_mrav = copy.deepcopy(
                    min(self.najbolji_mrav, najuspesniji_mrav_generacije))
            
            # pripremamo mrave za ponavljanje tako sto ih resetujemo
            for mrav in self.mravi:
                mrav.inicijalizuj_atribute()
                
    def stampaj_resenje(self):  
        if self.najbolji_mrav:
            try:
                # ako je instalirana matplot biblioteka koristi je za graficki prikaz
                import matplotlib.pyplot as plt 
                
                x, y = self.najbolji_mrav.koordinate()
                                
                # 1. tacku dodajemo ponovo da bi zatvorili konturu od poslednje ka 1. tacki
                x.append(x[0])
                y.append(y[0])
                
                plt.title('prikaz nadjenog resenja')
                plt.plot(x, y, 'g', marker='.', markerfacecolor='r', markersize=20, label='najbolje r.')
                plt.xlabel('x koordinata')
                plt.ylabel('y koordinata')
                plt.text(2, 95, 'najbolja duzina: %f.2' % self.najbolji_mrav.duzina_putanje)
                plt.legend()
                
                for grad in self.najbolji_mrav.putanja:
                    plt.text(grad.x, grad.y - 2.5, grad.ime)
                
                plt.axis([min(x) -1, max(x) + 1, min(y) -1, max(y) +1])             
             
                plt.show()
            except ImportError:
                # nema matplot biblioteke stampaj resenje rucno
                print "Resenje"
                print self.najbolji_mrav.putanja
                print "duzina puta je: ", self.najbolji_mrav.duzina_putanje
                
        
def main():
    try:
        import psyco # ako je instalirana biblioteka koristi je (ubrzava izvsenje koda)
        psyco.full()
    except ImportError:
        pass    
    
    gradovi = [Grad(1, 2, 'a'), Grad(29, 21, 'b'),Grad(100, 60, 'c'), Grad(12, 86, 'd'), 
            Grad(92, 46, 'e'), Grad(83, 38, 'f'), Grad(55, 36, 'g'), Grad(71, 99, 'h'),
            Grad(12, 41, 'i'), Grad(34, 48, 'j'), Grad(69, 33, 'k'), Grad(78, 10, 'l'),
            Grad(86, 68, 'm'), Grad(79, 27, 'n'), Grad(22, 69, 'o'), Grad(75, 55, 'p'),
            Grad(51, 68, 'r'), Grad(91, 23, 's'), Grad(22, 42, 't'), Grad(47, 80, 'u'), 
            Grad(60, 10, 'w'), Grad(91, 79, 'v'), Grad(5, 66, 'x'), Grad(42, 90, 'y'),
            Grad(23, 59, '1'), Grad(46, 83, '2'), Grad(93, 63, '3'), Grad(47, 17, '4'),
            Grad(53, 79, '5'), Grad(76, 23, '6'), Grad(91, 62, '7'), Grad(44, 97, '8'),
            ]
            
    simulator = MravljiSimulator(gradovi, ALFA=1.0, BETA = 2.0, RO=0.5)
    
    t0 = time.time()
    simulator.start(ponavljanja=50)
    print "Izvrsenje trajalo", time.time() - t0, "sekundi"
    
    simulator.stampaj_resenje()

if __name__ == "__main__":
    main()    
    

Gradovi koje sam ovde koristio su isti oni koje sam korisitio sa Simuliranim kaljenjem. Rastojanja gradova se unapred izračunaju da bi se program brže izvršavao. Brzo se dobijaju dobra rešenja. Otprilike da je putanja koju pronađe u proseku dugačka oko 530. A svaki 15-20 put ubode se i rešenje čiju sam sliku postavio. Za razliku od Simuliranog kaljena ovde ima mnoštvo parametara sa kojima se treba igrati da se nađe njihova dobra vrednost.