Stackelberg Bi-Level GA attempt

import random
import matplotlib
import numpy as np

matplotlib.use('TkAgg')
import matplotlib.pyplot as plt

from deap import base, algorithms
from deap import creator
from deap import tools

# konstanty genetického algorytmu
P_KRIZENI = 0.9 # pravdepodobnost krizeni
P_PROMENA = 0.1 # pravdepodobnost mutace
POPULATION_SIZE = 100 # maximalni pocet induvidui v populaci
MAX_GENERATIONS = 100 # maximalni pocet generaci

# konstanty programu
TYP_AUTA_DODAVKA = 'D'
TYP_AUTA_PIKUP = 'P'
SEZNAM_AUT = [TYP_AUTA_DODAVKA,TYP_AUTA_PIKUP]
CENA_DORUCENI = 80  # cena, kterou uctuje spolecnost zakaznikovi
ODMENA_DORUCENI_DODAVKOU = 30 # odmena, kterou ziskava kuryr na dodavce
ODMENA_DORUCENI_PIKAPEM = 45 # odmena, kterou ziskava kuryr na pikapu
NAKLAD_DORUCENI_DODAVKOU = 10 # palivo a mzda na jednu dodavku (kc/km, kc/h)
NAKLAD_DORUCENI_PIKAPEM = 15 # palivo a mzda na jeden pikap (kc/km, kc/h)
VMIN_DODAVKA, VMAX_DODAVKA = 1, 35 # minimalni a maximalni pocet baliku v dodavce
VMIN_PIKAP, VMAX_PIKAP = 1, 15 # minimalni a maximalni pocet baliku v pikapu
QMIN_DODAVKA, QMAX_DODAVKA = 0, 10 # minimalni a maximalni pocet dodavek u spolecnosti
QMIN_PIKAP, QMAX_PIKAP = 0, 15 # minimalni a maximalni pocet pikapu u spolecnosti

# registrace funkce obsahujici parametry, jez se budou optimalizovat pro auta
creator.create("FitnessMaxCars", base.Fitness, weights=(1.0,))
creator.create("IndividualCar", list, fitness=creator.FitnessMaxCars)

# registrace funkce obsahujici parametry, jez se budou optimalizovat pro spolecnosti
creator.create("FitnessMaxCompanies", base.Fitness, weights=(1.0,))
creator.create("IndividualCompany", list, fitness=creator.FitnessMaxCompanies)

# cilova funkce kuryra
def oneMaxFitnessKuryr(individual):
    zisk_kuryra, prinos_dodavky, prinos_pikapu = 0, 0, 0 # iniciace promennych
    typ_auta = individual[0] # urceni dle vstupu typu auta, na nemz jezdi
    pocet_doruceni = individual[1] # odvozeni dle vstupu poctu doruceni timto autem
    # urceni prinosu auta na zaklade jeho typu
    if typ_auta == TYP_AUTA_DODAVKA:
        prinos_dodavky = ODMENA_DORUCENI_DODAVKOU * pocet_doruceni
    elif typ_auta == TYP_AUTA_PIKUP:
        prinos_pikapu = ODMENA_DORUCENI_PIKAPEM * pocet_doruceni
    zisk_kuryra = prinos_dodavky + prinos_pikapu # secteni prinosu, (prinos + 0) je zisk kuryra
    return zisk_kuryra, # tuple

def oneMaxFitnessSpolecnost(individual):
    zisk_spolecnosti, prinos_dodavky, prinos_pikapu = 0, 0, 0 # iniciace promennych
    # Pocitame pocet aut kazdeho typ - potreba pro kontrolu omezeni
    count_dodavka = sum(1 for car in individual if car[0] == TYP_AUTA_DODAVKA)
    count_pikap = sum(1 for car in individual if car[0] == TYP_AUTA_PIKUP)

    # Je-li pocet aut dodavka nebo pikap vyssi nez pripustne cislo, vypisujeme pokutu na zisk spolecnosti
    if count_dodavka > QMAX_DODAVKA:
        zisk_spolecnosti -= 100000  # stanovme pokutu zisku
    if count_pikap > QMAX_PIKAP:
        zisk_spolecnosti -= 100000  # stanovme pokutu zisku

    for car in individual:
        typ_auta = car[0]
        pocet_doruceni = car[1]

        # pocitame zisk spolecnosti
        if typ_auta == TYP_AUTA_DODAVKA:
            prinos_dodavky = (CENA_DORUCENI - ODMENA_DORUCENI_DODAVKOU) * pocet_doruceni * NAKLAD_DORUCENI_DODAVKOU
        elif typ_auta == TYP_AUTA_PIKUP:
            prinos_pikapu = (CENA_DORUCENI - ODMENA_DORUCENI_PIKAPEM) * pocet_doruceni * NAKLAD_DORUCENI_PIKAPEM

        zisk_spolecnosti += prinos_dodavky + prinos_pikapu

    return zisk_spolecnosti,

# souprava knihovny DEAP k provedeni genetickeho algoritmu
toolbox = base.Toolbox()

# definujeme funkci pro tvorbu auta dle vzoru [typ, vytizeni]
def createIndividualCarTypeAndCargoVolume(type):
    if type == TYP_AUTA_DODAVKA:
        return [type, random.randint(VMIN_DODAVKA,VMAX_DODAVKA)]
    if type == TYP_AUTA_PIKUP:
        return [type, random.randint(VMIN_PIKAP,VMAX_PIKAP)]

# definujeme funkci pro tvorbu poctu aut spolecnosti dle vzoru [dodavky, pikapy]
def createIndividualCompanyCars():
    D = random.randint(QMIN_DODAVKA, QMAX_DODAVKA) # nahodny pocet aut typu Dodavka
    P = random.randint(QMIN_PIKAP, QMAX_PIKAP) # nahodny pocet aut typu Pikap
    dodavky = [createIndividualCarTypeAndCargoVolume(TYP_AUTA_DODAVKA) for _ in range(D)] # nahodne vytvarime vytizeni kazde dodavce
    pikapy = [createIndividualCarTypeAndCargoVolume(TYP_AUTA_PIKUP) for _ in range(P)] # nahodne vytvarime vytizeni kazdemu pikapu
    auta = dodavky + pikapy # sjednotime seznamy dodavek a pikapu na jeden seznam
    return auta

toolbox.register("carAssign", createIndividualCompanyCars) # registrujeme v DEAP funkci vytvareni aut vcetne nakladu, viz vyse
toolbox.register("individualCompanyCreator", tools.initIterate, creator.IndividualCompany, toolbox.carAssign) # registrujeme v DEAP funkci tvorby jednotlive spolecnosti
toolbox.register("populationCompanyCreator", tools.initRepeat, list, toolbox.individualCompanyCreator) # registrujeme funkci pro tvorbu populace spolecnosti

# toolbox.register("typeAndVolumeAssign", createIndividualCarTypeAndCargoVolume) # registrujeme funkci pro tvorbu typu auta a jeho nakladu
# toolbox.register("individualCarCreator", tools.initIterate, creator.IndividualCar, toolbox.typeAndVolumeAssign) # registrujeme funkci pro tvorbu jednotliveho auta
# toolbox.register("populationCarCreator", tools.initRepeat, list, toolbox.individualCarCreator) # registrujeme funci pro tvorbu populace aut

populationCompanies = toolbox.populationCompanyCreator(n=POPULATION_SIZE) # vytvarime populaci spolecnosti
# prochazime kazdou spolecnist a vytvarime populaci aut z nich
populationCars = [creator.IndividualCar(auto) for spolecnost in [toolbox.carAssign() for _ in range(POPULATION_SIZE)] for auto in spolecnost]
# populationCars = toolbox.populationCarCreator(n=POPULATION_SIZE) # vytvarime populaci aut

# funkce uzivatelskeho krizeni
def customCrossover(ind1, ind2):
    # krizeni pro typ auta, vymenime mista
    ind1[0], ind2[0] = ind2[0], ind1[0]

    # krizeni pro pocet baliku - bereme prumer
    ind1[1] = (ind1[1] + ind2[1]) // 2
    ind2[1] = ind1[1]
    # pred vracenim hodnot kontrolujeme na prekroceni poctu baliku omezeni dle typu auta, pouzijeme minimum z hodnoty nebo omezeni
    if ind1[0] == TYP_AUTA_PIKUP:
        ind1[1] = min(ind1[1], VMAX_PIKAP)
    else:
        ind2[1] = min(ind2[1], VMAX_DODAVKA)

    return ind1, ind2

# funkce uzivatelske mutace
def customMutation(individual, indpb):
    if random.random() < indpb: # pravdepodobnost mutace

        individual[0] = random.choice(SEZNAM_AUT) # mutace typu auta - nahodny vyber

    # mutace poctu doruceni (baliku)
    if random.random() < indpb:
        if individual[0] == TYP_AUTA_DODAVKA:
            individual[1] = random.randint(VMIN_DODAVKA, VMAX_DODAVKA)
        elif individual[0] == TYP_AUTA_PIKUP:
            individual[1] = random.randint(VMIN_PIKAP, VMAX_PIKAP)

    return individual,

# registrace operatoru genetickeho algoritmu pro kuryra
toolbox.register("evaluate", oneMaxFitnessKuryr)
toolbox.register("select", tools.selTournament, tournsize=3)
toolbox.register("mate", customCrossover) # uzivatelske krizeni aut
toolbox.register("mutate", customMutation, indpb=P_PROMENA) # uzivatelska mutace aut

# vypocteme prizpusobenost kazdeho individua v nase sformovane populaci kuryru
print("Ohodnoceni aut...")
fitnessValuesCars = list(map(oneMaxFitnessKuryr, populationCars)) # funkce vytvori seznam hodnot cilove funkce kuryru pro kazdeho individua v soucasne populaci
# vlastnosti 'values', ktera je ulozena ve tride FitnessMaxCars pro auta predame vypoctene hodnoty cilove funkce
for individual, fitnessValue in zip(populationCars, fitnessValuesCars):
    individual.fitness.values = fitnessValue

# registrujeme soupravu pro sledovani statistik
stats = tools.Statistics(key=lambda ind: ind.fitness.values)
stats.register("max", np.max)
stats.register("avg", np.mean)

# provadime geneticky algoritmus pro populaci aut
populationCars, logbook = algorithms.eaSimple(populationCars, toolbox,
                                          cxpb=P_KRIZENI,
                                          mutpb=P_PROMENA,
                                          ngen=MAX_GENERATIONS,
                                          stats=stats,
                                          verbose=True)

# vypocteme prizpusobenost kazdeho individua v nase populaci spolecnosti
print("Ohodnoceni spolecnosti...")
fitnessValuesCompanies = list(map(oneMaxFitnessSpolecnost, populationCompanies)) # funkce vytvori seznam hodnot cilove funkce spolecnosti pro kazdeho individua v soucasne populaci
# vlastnosti 'values', ktera je ulozena ve tride FitnessMaxCompanies pro spolecnosti predame vypoctene hodnoty cilove funkce
for individual, fitnessValue in zip(populationCompanies, fitnessValuesCompanies):
    individual.fitness.values = fitnessValue

# funkce uzivatelskeho krizeni
def customMutationCompany(individual, indpb):
    if random.random() < indpb: # pravdepodobnost mutace
        # nahodne vybereme auto
        car_index = random.randint(0, len(individual) - 1)
        car = individual[car_index]

        # overime na omezeni poctu aut dle typu pred mutaci
        count_dodavka = sum(1 for c in individual if c[0] == TYP_AUTA_DODAVKA)
        count_pikap = sum(1 for c in individual if c[0] == TYP_AUTA_PIKUP)

        # pokud pocet aut vyhovuje omezeni, provedeme mutaci
        if car[0] == TYP_AUTA_DODAVKA and count_dodavka < QMAX_DODAVKA:
            car[0] = random.choice(SEZNAM_AUT)
        elif car[0] == TYP_AUTA_PIKUP and count_pikap < QMAX_PIKAP:
            car[0] = random.choice(SEZNAM_AUT)

        # mutace poctu doruceni (baliku)
        if car[0] == TYP_AUTA_DODAVKA:
            car[1] = random.randint(VMIN_DODAVKA, VMAX_DODAVKA)
        elif car[0] == TYP_AUTA_PIKUP:
            car[1] = random.randint(VMIN_PIKAP, VMAX_PIKAP)

    return individual,

# odebrani operatoru genetickeho algoritmu pro kuryra
toolbox.unregister("evaluate")
toolbox.unregister("select")
toolbox.unregister("mate")
toolbox.unregister("mutate")

# registrace operatoru genetickeho algoritmu pro spolecnost
toolbox.register("evaluate", oneMaxFitnessSpolecnost)
toolbox.register("select", tools.selTournament, tournsize=3)
toolbox.register("mate", tools.cxOnePoint)
toolbox.register("mutate", customMutationCompany, indpb=P_PROMENA)

# provadime geneticky algoritmus pro spolecnost
populationCompanies, logbook2 = algorithms.eaSimple(populationCompanies, toolbox,
                                          cxpb=P_KRIZENI,
                                          mutpb=P_PROMENA,
                                          ngen=MAX_GENERATIONS,
                                          stats=stats,
                                          verbose=True)

# zjistime nejlepsiho individa (spolecnost) po evoluci
best_individual = tools.selBest(populationCompanies, 1)[0]
print("Nejlepsi spolecnost je ", best_individual)

# vypocteme pocet kazdeho typu aut v nejlepsi spolecnosti
count_dodavka = sum(1 for car in best_individual if car[0] == TYP_AUTA_DODAVKA)
count_pikap = sum(1 for car in best_individual if car[0] == TYP_AUTA_PIKUP)

print("Pocet dodavek = ", count_dodavka)
print("Pocet pikapu = ", count_pikap)

# urcime dva seznamu pro uchovavani statistiky prace algorytmu, pro auta
maxFitnessValues = logbook.select("max") # maximalni prizpusobeni v soucasne populaci
meanFitnessValues = logbook.select("avg") # prumerne prizpusobeni v soucasne populaci
# urcime dva seznamu pro uchovavani statistiky prace algorytmu, pro spolecnosti
maxFitnessValues2 = logbook2.select("max") # maximalni prizpusobeni v soucasne populaci
meanFitnessValues2 = logbook2.select("avg") # prumerne prizpusobeni v soucasne populaci

# auta
plot1 = plt.figure(1)
plt.plot(maxFitnessValues, color='grey')
plt.plot(meanFitnessValues, color='red')
plt.xlabel('Generace')
plt.ylabel('Maximální/průměrný zisk kurýra (přizpůsobivost)')
plt.title('Závislost maximální a střední přizpůsobivosti na generaci')

# spolecnosti
plot2 = plt.figure(2)
plt.plot(maxFitnessValues2, color='grey')
plt.plot(meanFitnessValues2, color='red')
plt.xlabel('Generace')
plt.ylabel('Maximální/průměrný zisk společnosti (přizpůsobivost)')
plt.title('Závislost maximální a střední přizpůsobivosti na generaci')
plt.show()