Genetic Algorithm Lab

# Binary GA for the ticket pricing problem (Chapter 7, Sections 7.3 and 7.4)
# Maximize profit(x) = -20x^2 + 6200x - 350000 for 75 <= x <= 235, with 8-bit chromosomes.
import numpy as np

profit = lambda x: -20 * x * x + 6200 * x - 350000
BITS, LB, UB = 8, 75, 235
to_bits = lambda x: [int(b) for b in format(x, f"0{BITS}b")]
to_int = lambda bits: int("".join(map(str, bits)), 2)
gray = lambda b: [b[0]] + [b[i - 1] ^ b[i] for i in range(1, len(b))]      # binary -> Gray (XOR)
def from_gray(g):                                                            # Gray -> binary
    b = [g[0]]
    for i in range(1, len(g)):
        b.append(b[-1] ^ g[i])
    return b

# Table 7.2 population and the roulette wheel of Tables 7.4 and 7.5
pop = [77, 203, 110, 145, 230]
f = np.array([profit(x) for x in pop], dtype=float)
p = f / f.sum()
print("selection probabilities:", np.round(p, 2), " cumulative:", np.round(np.cumsum(p), 2))

# linear ranking with selection pressure SP (rank 1 = worst, N = best)
SP, N = 1.5, len(pop)
rank = np.argsort(np.argsort(f)) + 1
p_rank = (2 - SP + 2 * (SP - 1) * (rank - 1) / (N - 1)) / N
print("rank-based probabilities:", np.round(p_rank, 3))

# crossover and mutation on parents 77 and 145 (Figures 7.13 to 7.16)
P1, P2 = to_bits(77), to_bits(145)
c1, c2 = P1[:3] + P2[3:], P2[:3] + P1[3:]
print("1-point crossover:", to_int(c1), to_int(c2), "profits", profit(to_int(c1)), profit(to_int(c2)))
mut = P1.copy()
for j in (2, 4, 6):
    mut[j] ^= 1
print("mutated 77 ->", to_int(mut), "profit", profit(to_int(mut)))


def bga(pop_size=20, n_gen=60, pc=0.7, pm=0.05, selection="tournament", elite=1,
        use_gray=False, seed=0):
    rng = np.random.default_rng(seed)
    encode = (lambda x: gray(to_bits(x))) if use_gray else to_bits
    decode = (lambda c: to_int(from_gray(c))) if use_gray else to_int
    fit = lambda c: profit(decode(c)) if LB <= decode(c) <= UB else -1e9   # infeasible: rejected
    P = [encode(int(x)) for x in rng.integers(LB, UB + 1, pop_size)]
    history = []
    for g in range(n_gen):
        F = np.array([fit(c) for c in P])
        order = np.argsort(-F)
        history.append(F[order[0]])

        def pick():
            if selection == "roulette":
                w = F - F.min() + 1e-9
                return P[rng.choice(pop_size, p=w / w.sum())]
            if selection == "rank":
                r = np.argsort(np.argsort(F)) + 1
                w = (2 - SP + 2 * (SP - 1) * (r - 1) / (pop_size - 1)) / pop_size
                return P[rng.choice(pop_size, p=w / w.sum())]
            if selection == "random":
                return P[rng.integers(pop_size)]
            k = rng.choice(pop_size, 3, replace=False)                    # tournament, k = 3
            return P[k[np.argmax(F[k])]]

        nxt = [P[i].copy() for i in order[:elite]]                          # elitism
        while len(nxt) < pop_size:
            a, b = pick().copy(), pick().copy()
            if rng.random() < pc:                                           # 1-point crossover
                cut = rng.integers(1, BITS)
                a, b = a[:cut] + b[cut:], b[:cut] + a[cut:]
            for child in (a, b):
                for j in range(BITS):                                       # bit-flip mutation
                    if rng.random() < pm:
                        child[j] ^= 1
                nxt.append(child)
        P = nxt[:pop_size]
    best = max(P, key=fit)
    return decode(best), fit(best), history


x, best_profit, hist = bga(seed=1)
print(f"optimal ticket price: ${x}, profit: ${best_profit:,.0f}")      # $155, $130,500

# Real-valued GA for curve fitting (Chapter 8, Section 8.2)
# Fit y = a x^3 + b x^2 + c x + d to five points generated by a, b, c, d = 2, -3, 4, 1.
import numpy as np

x = np.linspace(0, 5, 5)
y = 2 * x**3 - 3 * x**2 + 4 * x + 1
mse = lambda coeffs: np.mean((y - np.polyval(coeffs, x)) ** 2)
print("MSE of parent 1 [1 2 3 4]:", round(mse([1, 2, 3, 4]), 2))    # about 122.8
print("MSE of parent 2 [2 2 2 2]:", round(mse([2, 2, 2, 2]), 2))    # 3662.32


def single_arithmetic(p1, p2, k, alpha):
    c1, c2 = p1.copy(), p2.copy()
    c1[k] = alpha * p1[k] + (1 - alpha) * p2[k]
    c2[k] = (1 - alpha) * p1[k] + alpha * p2[k]
    return c1, c2


def simple_arithmetic(p1, p2, k, alpha):
    c1, c2 = p1.copy(), p2.copy()
    c1[k:] = alpha * p1[k:] + (1 - alpha) * p2[k:]
    c2[k:] = (1 - alpha) * p1[k:] + alpha * p2[k:]
    return c1, c2


def whole_arithmetic(p1, p2, alpha):
    return alpha * p1 + (1 - alpha) * p2, (1 - alpha) * p1 + alpha * p2


A = np.array([1.1, 4.3, 2.3, 6.7, 1.5, 0.6, 1.8, 2.5])
B = np.array([2.3, 3.2, 1.3, 4.2, 5.3, 2.2, 3.3, 4.2])
print("single (k = 6):", single_arithmetic(A, B, 5, 0.3)[0])        # gene 6 -> 1.72 (Figure 8.2)
print("whole:", np.round(whole_arithmetic(A, B, 0.3)[0], 2))         # 1.94 3.53 ... (Figure 8.4)
print("whole on [1 2 3 4], [2 2 2 2], alpha 0.2:", whole_arithmetic(np.array([1., 2, 3, 4]), np.array([2., 2, 2, 2]), 0.2)[0])


def real_ga(pop_size=60, n_gen=200, pm=0.3, sigma0=1.0, crossover="whole",
            mutation="gaussian", adaptive=True, lb=-10, ub=10, seed=0):
    rng = np.random.default_rng(seed)
    P = rng.uniform(lb, ub, (pop_size, 4))
    history = []
    for g in range(n_gen):
        F = np.array([mse(p) for p in P])
        history.append(F.min())
        sigma = sigma0 * (1 - 0.9 * g / n_gen) if adaptive else sigma0    # deterministic schedule
        elite = P[np.argsort(F)[:2]]
        children = [*elite]
        while len(children) < pop_size:
            i, j = (min(rng.choice(pop_size, 3, replace=False), key=lambda k: F[k]) for _ in range(2))
            alpha = rng.random()
            if crossover == "whole":
                c1, c2 = whole_arithmetic(P[i], P[j], alpha)
            elif crossover == "simple":
                c1, c2 = simple_arithmetic(P[i], P[j], rng.integers(4), alpha)
            else:
                c1, c2 = single_arithmetic(P[i], P[j], rng.integers(4), alpha)
            for c in (c1, c2):
                for d in range(4):
                    if rng.random() < pm:
                        if mutation == "gaussian":
                            c[d] += rng.normal(0, sigma)
                        elif mutation == "cauchy":
                            c[d] += sigma * rng.standard_cauchy()
                        else:                                                 # boundary (uniform)
                            c[d] = rng.uniform(lb, ub)
                children.append(np.clip(c, lb, ub))
        P = np.array(children[:pop_size])
    best = P[np.argmin([mse(p) for p in P])]
    return best, mse(best), history


best, err, hist = real_ga(seed=2)
print("a, b, c, d =", np.round(best, 3), " MSE =", round(err, 6))     # close to 2, -3, 4, 1

# Permutation-based GA (Chapter 8, Sections 8.3 and 8.6)
# Crossovers: PMX, edge (EC), order 1 (OX1), cycle (CX); mutations: insert, swap, inversion, scramble.
import numpy as np


def pmx(p1, p2, a, b):
    child = [None] * len(p1)
    child[a:b + 1] = p1[a:b + 1]
    for i in range(len(p1)):
        if child[i] is None:
            v = p2[i]
            while v in child:                     # follow the mapping until v is free
                v = p2[p1.index(v)]
            child[i] = v
    return child


def ox1(p1, p2, a, b):
    n, child = len(p1), [None] * len(p1)
    child[a:b + 1] = p1[a:b + 1]
    fill = [p2[(b + 1 + k) % n] for k in range(n) if p2[(b + 1 + k) % n] not in child]
    for k, pos in enumerate([(b + 1 + k) % n for k in range(n - (b - a + 1))]):
        child[pos] = fill[k]
    return child


def cx(p1, p2):
    n, c1, c2, seen, use_p1 = len(p1), [None] * len(p1), [None] * len(p1), set(), True
    for start in range(n):
        if start in seen:
            continue
        cycle, i = [], start
        while i not in cycle:                      # positions linked through the two parents
            cycle.append(i)
            i = p1.index(p2[i])
        seen.update(cycle)
        for i in cycle:
            c1[i], c2[i] = (p1[i], p2[i]) if use_p1 else (p2[i], p1[i])
        if len(cycle) > 1:                         # identical genes do not switch the parent
            use_p1 = not use_p1
    return c1, c2


def edge_crossover(p1, p2, start=1):
    n = len(p1)
    table = {c: [] for c in p1}
    for p in (p1, p2):
        for i, c in enumerate(p):
            table[c] += [p[i - 1], p[(i + 1) % n]]
    child, cur = [start], start
    while len(child) < n:
        for lst in table.values():                 # remove the current city everywhere
            while cur in lst:
                lst.remove(cur)
        cand = table[cur]
        common = sorted({c for c in cand if cand.count(c) > 1})
        if common:
            nxt = common[0]
        elif cand:                                 # shortest edge list, ties by smallest label
            nxt = min(set(cand), key=lambda c: (len(set(table[c])), c))
        else:
            nxt = min(c for c in p1 if c not in child)
        child.append(nxt)
        cur = nxt
    return child


P1, P2 = [1, 2, 3, 4, 5, 6, 7, 8], [1, 6, 7, 8, 5, 2, 4, 3]
print("PMX :", pmx(P1, P2, 2, 5))        # [1, 2, 3, 4, 5, 6, 8, 7]  (Figure 8.5)
print("EC  :", edge_crossover(P1, P2))   # [1, 3, 4, 2, 5, 6, 7, 8]  (Figure 8.7)
print("OX1 :", ox1(P1, P2, 2, 5))        # [8, 2, 3, 4, 5, 6, 1, 7]  (Figure 8.8)
print("CX  :", cx(P1, P2))               # [1, 2, 7, 8, 5, 6, 4, 3], [1, 6, 3, 4, 5, 2, 7, 8]
inversion = lambda t, i, j: t[:i] + t[i:j + 1][::-1] + t[j + 1:]
print("inversion:", inversion(P1, 1, 4))  # [1, 5, 4, 3, 2, 6, 7, 8]  (Figure 8.10)

# GA for the 20-city TSP (closed tour from New York City), OX1 + inversion mutation
cities = {
    "New York City": (40.72, -74.00), "Philadelphia": (39.95, -75.17), "Baltimore": (39.28, -76.62),
    "Charlotte": (35.23, -80.85), "Memphis": (35.12, -89.97), "Jacksonville": (30.32, -81.70),
    "Houston": (29.77, -95.38), "Austin": (30.27, -97.77), "San Antonio": (29.53, -98.47),
    "Fort Worth": (32.75, -97.33), "Dallas": (32.78, -96.80), "San Diego": (32.78, -117.15),
    "Los Angeles": (34.05, -118.25), "San Jose": (37.30, -121.87), "San Francisco": (37.78, -122.42),
    "Indianapolis": (39.78, -86.15), "Phoenix": (33.45, -112.07), "Columbus": (39.98, -82.98),
    "Chicago": (41.88, -87.63), "Detroit": (42.33, -83.05),
}
names = list(cities)
lat, lon = np.radians(np.array(list(cities.values()))).T
a = np.sin((lat[:, None] - lat) / 2) ** 2 + np.cos(lat[:, None]) * np.cos(lat) * np.sin((lon[:, None] - lon) / 2) ** 2
D = 2 * 6371.0088 * np.arcsin(np.sqrt(a))
n = len(names)
tour_len = lambda t: sum(D[t[k], t[(k + 1) % n]] for k in range(n))


def tsp_ga(pop_size=60, n_offspring=60, n_gen=500, seed=0):
    """(mu + lambda) GA: OX1 crossover, inversion mutation on every child, no duplicates."""
    rng = np.random.default_rng(seed)
    P = [[0] + list(rng.permutation(np.arange(1, n))) for _ in range(pop_size)]   # start at NYC
    F = [tour_len(t) for t in P]
    for g in range(n_gen):
        kids = []
        for _ in range(n_offspring):
            i, j = (min(rng.choice(pop_size, 2, replace=False), key=lambda k: F[k]) for _ in range(2))
            a_, b_ = sorted(rng.choice(np.arange(1, n), 2, replace=False))
            child = [0] + ox1(P[i][1:], P[j][1:], a_ - 1, b_ - 1)                  # keep NYC first
            u, v = sorted(rng.choice(np.arange(1, n), 2, replace=False))
            kids.append(inversion(child, u, v))
        pool = {tuple(t): f for t, f in zip(P + kids, F + [tour_len(t) for t in kids])}
        best = sorted(pool.items(), key=lambda kv: kv[1])[:pop_size]          # survivor selection
        P, F = [list(t) for t, _ in best], [f for _, f in best]
    return P[0], F[0]


best, L = tsp_ga(seed=3)          # some seeds stop at 11,005.2 km, a near-optimal tour
print(" -> ".join(names[i] for i in best))
print(f"Tour length: {L:.1f} km")   # optimum 10,934.8 km

# NSGA-II from scratch for the manufacturing problem (Chapter 8, Section 8.4)
# x1, x2: units of P1 on M1, M2; y1, y2: units of P2 on M1, M2.
# minimize f1 = total cost and f2 = |time on M1 - time on M2|, s.t. x1 + x2 >= 80, y1 + y2 >= 300.
import numpy as np
import matplotlib.pyplot as plt

LB, UB = np.zeros(4), np.array([100, 100, 500, 500])


def evaluate(v):
    x1, x2, y1, y2 = v
    f = np.array([100 * x1 + 150 * x2 + 200 * y1 + 50 * y2, abs((2 * x1 + 4 * y1) - (3 * x2 + y2))])
    cv = max(0, 80 - x1 - x2) + max(0, 300 - y1 - y2)          # total constraint violation
    return f, cv


def dominates(a, b):                                           # constraint-domination (Deb)
    (fa, ca), (fb, cb) = a, b
    if ca != cb:
        return ca < cb
    return np.all(fa <= fb) and np.any(fa < fb)


def nondominated_sort(E):
    n = len(E)
    S, cnt, fronts = [[] for _ in range(n)], np.zeros(n, int), [[]]
    for p in range(n):
        for q in range(n):
            if dominates(E[p], E[q]):
                S[p].append(q)
            elif dominates(E[q], E[p]):
                cnt[p] += 1
        if cnt[p] == 0:
            fronts[0].append(p)
    while fronts[-1]:
        nxt = []
        for p in fronts[-1]:
            for q in S[p]:
                cnt[q] -= 1
                if cnt[q] == 0:
                    nxt.append(q)
        fronts.append(nxt)
    return fronts[:-1]


def crowding(F):
    d = np.zeros(len(F))
    for m in range(F.shape[1]):
        o = np.argsort(F[:, m]); d[o[0]] = d[o[-1]] = np.inf
        span = F[o[-1], m] - F[o[0], m] or 1
        d[o[1:-1]] += (F[o[2:], m] - F[o[:-2], m]) / span
    return d


def sbx(a, b, rng, eta=15):
    u = rng.random(4)
    beta = np.where(u <= .5, (2 * u) ** (1 / (eta + 1)), (1 / (2 * (1 - u))) ** (1 / (eta + 1)))
    return .5 * ((1 + beta) * a + (1 - beta) * b), .5 * ((1 - beta) * a + (1 + beta) * b)


def poly_mutation(v, rng, eta=20, pm=0.25):
    v = v.copy()
    for k in range(4):
        if rng.random() < pm:
            u = rng.random()
            delta = (2 * u) ** (1 / (eta + 1)) - 1 if u < .5 else 1 - (2 * (1 - u)) ** (1 / (eta + 1))
            v[k] += delta * (UB[k] - LB[k])
    return np.clip(v, LB, UB)


def nsga2(pop_size=40, n_gen=200, seed=1):
    rng = np.random.default_rng(seed)
    P = rng.uniform(LB, UB, (pop_size, 4))
    for g in range(n_gen):
        E = [evaluate(v) for v in P]
        rank = np.zeros(pop_size, int)
        for r, fr in enumerate(nondominated_sort(E)):
            rank[fr] = r
        cd = crowding(np.array([e[0] for e in E]))
        def tour():                                            # binary tournament: rank, then crowding
            i, j = rng.choice(pop_size, 2, replace=False)
            return P[i] if (rank[i], -cd[i]) < (rank[j], -cd[j]) else P[j]
        Q = []
        while len(Q) < pop_size:
            c1, c2 = sbx(tour(), tour(), rng)
            Q += [poly_mutation(c1, rng), poly_mutation(c2, rng)]
        R = np.vstack([P, Q]); ER = [evaluate(v) for v in R]
        new = []
        for fr in nondominated_sort(ER):                       # elitist survivor selection
            if len(new) + len(fr) <= pop_size:
                new += fr
            else:
                cdf = crowding(np.array([ER[i][0] for i in fr]))
                new += [fr[i] for i in np.argsort(-cdf)[:pop_size - len(new)]]
                break
        P = R[new]
    E = [evaluate(v) for v in P]
    front = nondominated_sort(E)[0]
    return P[front], np.array([E[i][0] for i in front])


X, F = nsga2()
o = np.argsort(F[:, 0]); X, F = X[o], F[o]
nF = (F - F.min(0)) / (F.max(0) - F.min(0) + 1e-12)            # normalize with ideal and nadir points
w = np.array([0.2, 0.8])
i = np.argmin(np.max(nF / (1 / w), axis=1))                   # augmented scalarization (ASF) pick
print("best trade-off: x1, x2, y1, y2 =", np.round(X[i]).astype(int), " cost =", round(F[i, 0], 1), " time difference =", round(F[i, 1], 1))
plt.scatter(F[:, 0], F[:, 1], facecolors="none", edgecolors="b"); plt.scatter(*F[i], marker="x", c="r", s=200)
plt.xlabel("total production cost"); plt.ylabel("difference in production times"); plt.show()

# PID tuning with a real-valued GA (Chapter 8, Section 8.7)
# Closed loop: T(s) = (Kd s^2 + Kp s + Ki) / (s^3 + (Kd + 10) s^2 + (Kp + 20) s + Ki)
import numpy as np
import matplotlib.pyplot as plt


def step_response(Kp, Ki, Kd, t_end=4.0, dt=0.002):
    """Simulate the plant 1/(s^2 + 10 s + 20) under PID control for a unit step."""
    n = int(t_end / dt)
    y = v = integ = 0.0
    e_prev = 1.0
    t, out = np.arange(n) * dt, np.zeros(n)
    for k in range(n):
        e = 1.0 - y
        integ += e * dt
        u = Kp * e + Ki * integ + Kd * (e - e_prev) / dt if k else Kp * e + Ki * integ
        e_prev = e
        a = u - 10 * v - 20 * y                   # y'' + 10 y' + 20 y = u
        v += a * dt; y += v * dt
        out[k] = y
    return t, out


def itae(params):
    t, y = step_response(*params)
    if not np.all(np.isfinite(y)) or np.abs(y).max() > 1e3:
        return 1e9
    return np.trapezoid(t * np.abs(1 - y), t)


def tune_pid(pop_size=30, n_gen=40, seed=0):
    rng = np.random.default_rng(seed)
    P = rng.uniform(0, 100, (pop_size, 3))
    for g in range(n_gen):
        F = np.array([itae(p) for p in P])
        sigma = 20 * (1 - 0.9 * g / n_gen)
        nxt = [P[np.argmin(F)]]
        while len(nxt) < pop_size:
            i, j = (min(rng.choice(pop_size, 3, replace=False), key=lambda k: F[k]) for _ in range(2))
            alpha = rng.random()
            child = alpha * P[i] + (1 - alpha) * P[j]                 # whole arithmetic crossover
            child += rng.normal(0, sigma, 3) * (rng.random(3) < 0.3)   # Gaussian mutation
            nxt.append(np.clip(child, 0, 100))
        P = np.array(nxt)
    best = min(P, key=itae)
    return best, itae(best)


(Kp, Ki, Kd), J = tune_pid()
t, y = step_response(Kp, Ki, Kd)
rise = t[np.argmax(y >= .9)] - t[np.argmax(y >= .1)]
print(f"Kp = {Kp:.1f}, Ki = {Ki:.1f}, Kd = {Kd:.1f}, ITAE = {J:.4f}, overshoot = {max(0, y.max() - 1) * 100:.1f}%, rise time = {rise:.3f} s")
plt.plot(t, y); plt.axhline(1, ls="--", c="gray"); plt.xlabel("time (s)"); plt.title("Step response"); plt.show()