import matplotlib.pyplot as plt import numpy as np import random import pandas as pd

遺傳算法求函數(shù)最值

遺傳算法的特點(diǎn)有無(wú)須可導(dǎo)，可優(yōu)化離散異構(gòu)數(shù)據(jù)。

參考：

• 遺傳算法python實(shí)現(xiàn)
• 科學(xué)前沿：Genetic Algorithm - 遺傳算法究竟怎么回事

問(wèn)題定義

如下圖所示，

在區(qū)間[-1, 2]上有很多極大值和極小值，如果要求其最大值或最小值，很多單點(diǎn)優(yōu)化的方法（梯度下降等）就不適合，這種情況下就可以用遺傳算法（Genetic Algorithm）。def fun(x):return x * np.sin(10*np.pi*x) + 2Xs = np.linspace(-1, 2, 100) plt.plot(Xs, fun(Xs))

初始化原始種群

遺傳算法的一個(gè)特點(diǎn)是同時(shí)優(yōu)化一批解（種群），例如下面在[-1, 2]范圍內(nèi)隨機(jī)生成

個(gè)點(diǎn)：np.random.seed(0)# 初始化原始種群 def ori_popular(num, _min=-1, _max=2):return np.random.uniform(_min, _max, num) # 范圍[-1, 2)#__TEST__ population = ori_popular(10)for pop, fit in zip(population, fun(population)):print("x=%5.2f, fit=%.2f"%(pop, fit))plt.plot(Xs, fun(Xs)) plt.plot(population, fun(population), '*') plt.show()>>> x= 0.65, fit=2.64 x= 1.15, fit=0.87 x= 0.81, fit=2.21 x= 0.63, fit=2.56 x= 0.27, fit=2.21 x= 0.94, fit=1.13 x= 0.31, fit=1.88 x= 1.68, fit=3.17 x= 1.89, fit=2.53 x= 0.15, fit=1.85

上圖顯示我們隨機(jī)選的

個(gè)點(diǎn)離最大值（3.8左右）差距還挺遠(yuǎn)，下面用GA算法看能否求到最優(yōu)解。

編碼

編碼，也就是由表現(xiàn)型到基因型，性征到染色體。

二進(jìn)制編碼的缺點(diǎn)：對(duì)于一些連續(xù)函數(shù)的優(yōu)化問(wèn)題，由于其隨機(jī)性使得其局部搜索能力較差，如對(duì)于一些高精度的問(wèn)題，當(dāng)解迫近于最優(yōu)解后，由于其變異后表現(xiàn)型變化很大，不連續(xù)，所以會(huì)遠(yuǎn)離最優(yōu)解，達(dá)不到穩(wěn)定。而格雷碼能有效地防止這類現(xiàn)象。

TODO:

• [ ] 用2**18擴(kuò)展似乎會(huì)損失精度，準(zhǔn)備用10000代替

下面分別將1, 10, 0轉(zhuǎn)成二進(jìn)制編碼，注意浮點(diǎn)數(shù)0.1無(wú)法轉(zhuǎn)換：

print("1:", bin(1)) print("10:", bin(10)) print("0:", bin(0))try:print("0.1:", bin(0.1)) except Exception as E:print("Exception: {}".format(type(E).__name__), E)>>> 1: 0b1 10: 0b1010 0: 0b0 Exception: TypeError 'float' object cannot be interpreted as an integer

為了能編碼浮點(diǎn)數(shù)，需要擴(kuò)大倍數(shù)轉(zhuǎn)成整數(shù)：

X = [1, 0.1, 0.002]print("2**18 =", 2**18, 'n') for x in X:tx = int(x * 2**18)print("%.4f => %6d => %s"%(x, tx, bin(tx)))2**18 = 262144 >>> 1.0000 => 262144 => 0b1000000000000000000 0.1000 => 26214 => 0b110011001100110 0.0020 => 524 => 0b1000001100for x in X:tx = int(x * 2**18)ex = bin(tx)dx = int(ex, 2) / 2**18print("%25s => %6d => %.4f"%(ex, tx, dx))>>> 0b1000000000000000000 => 262144 => 1.00000b110011001100110 => 26214 => 0.10000b1000001100 => 524 => 0.0020

需要注意的是，上面用%.4f進(jìn)行了四舍五入，實(shí)際上用2**18=262144來(lái)放大編碼會(huì)有精度的損失，例如：

x = 0.1 tx = int(x * 2**18) ex = bin(tx) int('0b110011001100110', 2) / (2**18)>>> 0.09999847412109375def encode(popular, _min=-1, _max=2, scale=2**18, bin_len=18): # popular應(yīng)該是float類型的列表"""bin_len: 染色體序列長(zhǎng)度"""norm_data = (popular-_min) / (_max-_min) * scalebin_data = np.array([np.binary_repr(x, width=bin_len) for x in norm_data.astype(int)]) # 轉(zhuǎn)成二進(jìn)制編碼return bin_datachroms = encode(population) for pop, chrom, fit in zip(population, chroms, fun(population)):print("x=%.2f, chrom=%s, fit=%.2f"%(pop, chrom, fit))>>> x=0.65, chrom=100011000111111100, fit=2.64 x=1.15, chrom=101101110001011010, fit=0.87 x=0.81, chrom=100110100100111010, fit=2.21 x=0.63, chrom=100010110111110101, fit=2.56 x=0.27, chrom=011011000111010010, fit=2.21 x=0.94, chrom=101001010101100101, fit=1.13 x=0.31, chrom=011100000000010110, fit=1.88 x=1.68, chrom=111001000100101100, fit=3.17 x=1.89, chrom=111101101011001010, fit=2.53 x=0.15, chrom=011000100010100100, fit=1.85

解碼和適應(yīng)度函數(shù)

通過(guò)基因（染色體）得到個(gè)體的適應(yīng)度值，也就是評(píng)判當(dāng)前種群每個(gè)個(gè)體（解）表現(xiàn)好壞，要對(duì)編碼解碼后才能代入函數(shù)。

注意，因?yàn)楹瘮?shù)

的結(jié)果就是大小，所以直接用來(lái)當(dāng)適應(yīng)度函數(shù)。def decode(popular_gene, _min=-1, _max=2, scale=2**18):return np.array([(np.int(x, base=2)/scale*3)+_min for x in popular_gene])fitness = fun(decode(chroms))for pop, chrom, dechrom, fit in zip(population, chroms, decode(chroms), fitness):print("x=%5.2f, chrom=%s, dechrom=%.2f, fit=%.2f"%(pop, chrom, dechrom, fit))>>> x= 0.65, chrom=100011000111111100, dechrom=0.65, fit=2.64 x= 1.15, chrom=101101110001011010, dechrom=1.15, fit=0.87 x= 0.81, chrom=100110100100111010, dechrom=0.81, fit=2.21 x= 0.63, chrom=100010110111110101, dechrom=0.63, fit=2.56 x= 0.27, chrom=011011000111010010, dechrom=0.27, fit=2.21 x= 0.94, chrom=101001010101100101, dechrom=0.94, fit=1.13 x= 0.31, chrom=011100000000010110, dechrom=0.31, fit=1.88 x= 1.68, chrom=111001000100101100, dechrom=1.68, fit=3.17 x= 1.89, chrom=111101101011001010, dechrom=1.89, fit=2.53 x= 0.15, chrom=011000100010100100, dechrom=0.15, fit=1.85

這里將最小的適應(yīng)度調(diào)整到0.000001，主要為了防止負(fù)數(shù)。

fitness = fitness - fitness.min() + 0.000001 # 保證所有的都為正

選擇與交叉

選擇就是淘汰不好的染色體（解），選擇方式有很多，這里用輪牌賭。

交叉就是讓染色體相互交換基因，方式有很多，這里在染色體中間進(jìn)行交叉，交叉概率默認(rèn)為0.66。

def SelectandCrossover(chroms, fitness, prob=0.6):probs = fitness/np.sum(fitness) # 各個(gè)個(gè)體被選擇的概率probs_cum = np.cumsum(probs) # 概率累加分布each_rand = np.random.uniform(size=len(fitness))#print(np.round(fitness, 2), np.round(probs, 2), np.round(probs_cum, 2), sep='n')#print([np.where(probs_cum > rand)[0][0] for rand in each_rand])# 根據(jù)隨機(jī)概率選擇出新的基因編碼newX = np.array([chroms[np.where(probs_cum > rand)[0][0]] for rand in each_rand])# 繁殖，隨機(jī)配對(duì)（概率為0.6)pairs = np.random.permutation(int(len(newX)*prob//2*2)).reshape(-1, 2)center = len(newX[0])//2for i, j in pairs:# 在中間位置交叉x, y = newX[i], newX[j]newX[i] = x[:center] + y[center:]newX[j] = y[:center] + x[center:]#print(x, y, newX[i], '', sep='n')return newX chroms = SelectandCrossover(chroms, fitness)dechroms = decode(chroms) fitness = fun(dechroms)for gene, dec, fit in zip(chroms, dechroms, fitness):print("chrom=%s, dec=%5.2f, fit=%.2f"%(gene, dec, fit))fig, (axs1, axs2) = plt.subplots(1, 2, figsize=(14, 5)) axs1.plot(Xs, fun(Xs)) axs1.plot(population, fun(population), 'o') axs2.plot(Xs, fun(Xs)) axs2.plot(dechroms, fitness, '*') plt.show()>>> chrom=111101101000010110, dec= 1.89, fit=2.64 chrom=011100000011001010, dec= 0.31, fit=1.86 chrom=011100000010100100, dec= 0.31, fit=1.86 chrom=011000100000010110, dec= 0.15, fit=1.85 chrom=100011000111111100, dec= 0.65, fit=2.64 chrom=100011000111111100, dec= 0.65, fit=2.64 chrom=100011000111111100, dec= 0.65, fit=2.64 chrom=111101101011001010, dec= 1.89, fit=2.53 chrom=111001000100101100, dec= 1.68, fit=3.17 chrom=111101101011001010, dec= 1.89, fit=2.53

變異

如果現(xiàn)在我們從生物學(xué)的角度來(lái)看這個(gè)問(wèn)題，那么請(qǐng)問(wèn)：由上述過(guò)程產(chǎn)生的后代是否有和其父母一樣的性狀呢？答案是否。在后代的生長(zhǎng)過(guò)程中，它們體內(nèi)的基因會(huì)發(fā)生一些變化，使得它們與父母不同。這個(gè)過(guò)程我們稱為「變異」，它可以被定義為染色體上發(fā)生的隨機(jī)變化，正是因?yàn)樽儺?#xff0c;種群中才會(huì)存在多樣性。

def Mutate(chroms:np.array):prob = 0.1clen = len(chroms[0])m = {'0':'1', '1':'0'}newchroms = []each_prob = np.random.uniform(size=len(chroms))for i, chrom in enumerate(chroms):if each_prob[i] < prob:pos = np.random.randint(clen) chrom = chrom[:pos] + m[chrom[pos]] + chrom[pos+1:]newchroms.append(chrom)return np.array(newchroms)newchroms = Mutate(chroms)def PltTwoChroms(chroms1, chroms2, fitfun):Xs = np.linspace(-1, 2, 100)fig, (axs1, axs2) = plt.subplots(1, 2, figsize=(14, 5))dechroms = decode(chroms1)fitness = fitfun(dechroms)axs1.plot(Xs, fitfun(Xs))axs1.plot(dechroms, fitness, 'o')dechroms = decode(chroms2)fitness = fitfun(dechroms)axs2.plot(Xs, fitfun(Xs))axs2.plot(dechroms, fitness, '*')plt.show()PltTwoChroms(chroms, newchroms, fun)

完整運(yùn)行

上面我們分解了每一步的運(yùn)算，下面我們用GA算法完整地迭代求出最優(yōu)解，這里種群個(gè)數(shù)為100，迭代次數(shù)為1000：

np.random.seed(0) population = ori_popular(100) chroms = encode(population)for i in range(1000):fitness = fun(decode(chroms))fitness = fitness - fitness.min() + 0.000001 # 保證所有的都為正newchroms = Mutate(SelectandCrossover(chroms, fitness))if i % 300 == 1:PltTwoChroms(chroms, newchroms, fun)chroms = newchromsPltTwoChroms(chroms, newchroms, fun)

深入思考

其實(shí)GA每一步都有各種各樣的算法，包括超參數(shù)（種群個(gè)數(shù)、交叉概率、變異概率），沒(méi)有一套組合適合所有的情況，需要根據(jù)實(shí)際情況來(lái)調(diào)整。例如我們?cè)谶M(jìn)行“選擇與交叉”時(shí)，可以考慮保留5%最優(yōu)秀的（因?yàn)橛幸欢ǜ怕时惶蕴?#xff09;，甚至這5%都不去交叉。

背包問(wèn)題

背包問(wèn)題是一個(gè)典型的優(yōu)化組合問(wèn)題：有一個(gè)容量固定的背包，已知

個(gè)物品的體積及其價(jià)值，我們?cè)撨x擇哪些物品裝入背包使得在其容量約束限制之內(nèi)所裝物品的價(jià)值最大？比如：背包容量限制為，有個(gè)物品的體積為，其價(jià)值則分別是。

背包問(wèn)題在計(jì)算理論中屬于NP-完全問(wèn)題，傳統(tǒng)的算法是動(dòng)態(tài)規(guī)劃法，可求出最優(yōu)總價(jià)值為

，其時(shí)間復(fù)雜度對(duì)于輸入是指數(shù)級(jí)的。而如果用貪心算法，計(jì)算物品性價(jià)比依次為，依次選擇物品和，總價(jià)值只能達(dá)到，這并非最優(yōu)解。

現(xiàn)在我們?cè)囉眠z傳算法求解這個(gè)問(wèn)題。

編碼

將待求解的各量表示成長(zhǎng)度為n的二進(jìn)制串。例如：

代表一個(gè)解，它表示將第、號(hào)物體放入背包中，其它的物體則不放入。

生成初始種群

本例中，群體規(guī)模的大小取為

，即群體由個(gè)個(gè)體組成，每個(gè)個(gè)體可通過(guò)隨機(jī)的方法產(chǎn)生。例如：np.random.seed(0) X = np.random.randint(0, 2, size=(4,4)) X>>> array([[0, 1, 1, 0],[1, 1, 1, 1],[1, 1, 1, 0],[0, 1, 0, 0]])

種群適應(yīng)度計(jì)算

對(duì)每個(gè)個(gè)體，計(jì)算所有物體的體積totalSize，所有物體的價(jià)值totalValue，以及個(gè)體適應(yīng)度f(wàn)it。若所有物體的總體積小于背包容量C，fit就是totalValue；否則應(yīng)進(jìn)行懲罰，fit應(yīng)減去alpha?(totalSize?C)，這里我們?nèi)lpha=2，得到如下結(jié)果：

each_size = np.array([2, 3, 4, 5]) each_value = np.array([3, 4, 5, 7]) C = 9 alpha = 2each_fit = np.zeros_like(each_size) for i, x in enumerate(X):total_size = (x*each_size).sum()total_value = (x*each_value).sum()fit = total_value if (total_size <= C) else total_value - alpha*(total_size - C)each_fit[i] = fitprint("x=%s, total_size=%d, total_value=%d, fit=%d"%(x, total_size, total_value, fit))>>> x=[0 1 1 0], total_size=7, total_value=9, fit=9 x=[1 1 1 1], total_size=14, total_value=19, fit=9 x=[1 1 1 0], total_size=9, total_value=12, fit=12 x=[0 1 0 0], total_size=3, total_value=4, fit=4

選擇

根據(jù)染色體適應(yīng)度，我們可以采用輪盤賭進(jìn)行染色體的選擇，即與適應(yīng)度成正比的概率來(lái)確定每個(gè)個(gè)體復(fù)制到下一代群體中的數(shù)量。具體操作過(guò)程如下：

先計(jì)算出群體中所有個(gè)體的適應(yīng)度的總和。

然后用每個(gè)個(gè)體適應(yīng)度除以總和，結(jié)果就是個(gè)體被遺傳到下一代群體中的概率。

每個(gè)概率值組成一個(gè)區(qū)間，全部概率值之和為。

最后再產(chǎn)生一個(gè)[0, 1]之間的隨機(jī)數(shù)，依據(jù)該隨機(jī)數(shù)出現(xiàn)在上述哪一個(gè)區(qū)間內(nèi)來(lái)確定各個(gè)個(gè)體被選中的次數(shù)。

fit_sum = each_fit.sum() probs = each_fit / fit_sum probs_cum = probs.cumsum()for i, x in enumerate(X):print('probility of %s is %.2f%%'%(x, probs[i]*100)) print('cumulative sum of the probilities:', np.round(probs_cum, 4)*100)>>> probility of [0 1 1 0] is 26.47% probility of [1 1 1 1] is 26.47% probility of [1 1 1 0] is 35.29% probility of [0 1 0 0] is 11.76% cumulative sum of the probilities: [ 26.47 52.94 88.24 100. ]new_num = 4 each_rand = np.random.uniform(size=4)np.round(each_rand, 2)>>> array([0.96, 0.38, 0.79, 0.53])

我們隨機(jī)產(chǎn)生了

個(gè)[0-1)之間的數(shù)each_rand，用于在probs_cum中符合的區(qū)間上選中對(duì)應(yīng)的基因編碼，分別選中了第88.24~100、26.47~52.94、52.94~88.24和26.47~52.94上的區(qū)間，所以新的基因編碼為：newX = np.array([ X[np.where(probs_cum > rand)[0][0]] for rand in each_rand])newX>>> array([[0, 1, 0, 0],[1, 1, 1, 1],[1, 1, 1, 0],[1, 1, 1, 1]])

沒(méi)想到概率最大的[1 1 1 0]反而被淘汰了。

交叉

采用單點(diǎn)交叉的方法，先對(duì)群體進(jìn)行隨機(jī)配對(duì)，其次隨機(jī)設(shè)置交叉點(diǎn)位置，最后再相互交換配對(duì)染色體之間的部分基因。

突變

我們采用基本位變異的方法來(lái)進(jìn)行變異運(yùn)算，首先確定出各個(gè)個(gè)體的基因變異位置，然后依照某一概率將變異點(diǎn)的原有基因值取反。若突變之后最大的適應(yīng)值個(gè)體已經(jīng)符合期望的價(jià)值則停止，否則返回第3步繼續(xù)迭代。這里我們?cè)O(shè)突變的概率為

：mutation_prob = 0.1 for i, x in enumerate(newX):if np.random.uniform() < mutation_prob:pos = np.random.randint(0, 4)newX[i][pos] = abs(x[pos] - 1) newX>>> array([[0, 1, 0, 0],[1, 1, 1, 1],[1, 0, 1, 0],[1, 1, 1, 1]])

總結(jié)

以上是生活随笔為你收集整理的python 求函数最大值_遗传算法与Python图解的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。