經典遺傳算法(SGA)解01揹包問題的原理及其python代碼實現

1.揹包問題

揹包問題（knapsack problem）是指從多種物品（項目）中選擇幾件物品轉滿揹包。假設存在n個不同的物品，對於物品j，其重量爲 $w_j$wj​，價值爲 $c_j$cj​，W是揹包承受的最大重量，揹包問題就是要在不超過揹包承受重量的前提下，使裝入揹包的物品的價值最大。

1.1簡單約束的揹包問題

揹包問題是理論上的NP-Hard問題，目前還沒有可求最優解的多項式時間算法。但很多情況下，採用遺傳算法在短時間內可以求的較好的解，定義 $x_j$xj​爲0-1決策變量，如果項目j被選擇，則 $x_j=1$xj​=1，否則 $x_j=0$xj​=0。下面描述一個簡單的帶有約束條件的揹包問題模型。
編號：
j：物品編號（j=1,2,…,n）
參數：
n：物品個數
W：揹包承受的最大重量
$w_j$wj​：物品j的重量
$c_j$cj​：物品j的價值
決策變量：
$x_j$xj​：物品j的決策變量，若物品j被選擇放入揹包，則 $x_j=1$xj​=1，否則 $x_j=0$xj​=0。
數學模型爲：
$s.t. \begin{matrix} max f(x)=\sum_{j=1}^{n}c_jx_j \\g(x)=\sum_{j=1}^{n}w_jx_j\leqslant W \\x_j=0/1,j=1,2,...,n \end{matrix}$s.t.maxf(x)=∑j=1n​cj​xj​g(x)=∑j=1n​wj​xj​⩽Wxj​=0/1,j=1,2,...,n​

1.2多選擇的揹包問題

該問題描述爲在若干類物品中有約束地選擇不同物品放入揹包使總價值最大，其數學模型如下：
編號：
i：類的編號（i=1,2,…,m）
j：物品編號（j=1,2,…, $n_i$ni​）
參數：
m：類的數量
$n_i$ni​：第i類中物品的數量
W：揹包承受的最大重量
$w_{ij}$wij​：第i類中第j個物品的重量
$c_{ij}$cij​：第i類中第j個物品的價值
決策變量：
$x_{ij}$xij​：第i類中第j個物品的決策變量，若第i類中第j個物品被選擇放入揹包，則 $x_{ij}=1$xij​=1，否則 $x_{ij}=0$xij​=0。
數學模型爲：
$s.t. \begin{matrix} max f(x)=\sum_{i=1}^{m}\sum_{j=1}^{n_i}c_{ij}x_{ij} \\g_w(x)=\sum_{i=1}^{m}\sum_{j=1}^{n_i}w_{ij}x_{ij}\leqslant W \\g_i(x)=\sum_{j=1}^{n_i}x_{ij}=1,\forall i \\x_j=0/1,\forall i,j \end{matrix}$s.t.maxf(x)=∑i=1m​∑j=1ni​​cij​xij​gw​(x)=∑i=1m​∑j=1ni​​wij​xij​⩽Wgi​(x)=∑j=1ni​​xij​=1,∀ixj​=0/1,∀i,j​

1.3多約束的揹包問題

這個問題是從多個類別中選擇物品，放入多個揹包中，在確保滿足揹包承重條件的前提下，使放入揹包的物品的總價值最大。表達式如下所示。
編號：
i：類的編號（i=1,2,…,m）
j：物品編號（j=1,2,…, $n_i$ni​）
參數：
m：類的數量
$n_i$ni​：第i類中物品的數量
$W_i$Wi​：揹包i承受的最大重量
$w_{ij}$wij​：第i類中第j個物品的重量
$c_{ij}$cij​：第i類中第j個物品的價值
決策變量：
$x_{ij}$xij​：第i類中第j個物品的決策變量，若第i類中第j個物品被選擇放入揹包i，則 $x_{ij}=1$xij​=1，否則 $x_{ij}=0$xij​=0。
數學模型爲：
$s.t. \begin{matrix} max f(x)=\sum_{i=1}^{m}\sum_{j=1}^{n_i}c_{ij}x_{ij} \\\sum_{j=1}^{n_i}w_{ij}x_{ij}\leqslant W_i,i=1,2,...,m \\x_j=0/1,\forall i,j \end{matrix}$s.t.maxf(x)=∑i=1m​∑j=1ni​​cij​xij​∑j=1ni​​wij​xij​⩽Wi​,i=1,2,...,mxj​=0/1,∀i,j​
本次實驗我們解決1.1問題，其它的問題讀者借鑑1.1問題的解決思路即可。
下面我們根據此問題介紹經典的遺傳算法。

2. 經典遺傳算法

2.1 遺傳算法概述

遺傳算法是一種基於生物遺傳機理，即生物進化（自然淘汰，交叉，變異等）現象的隨機搜索算法，它通過計算機模擬生物進化過程來進行搜索和進化，最後尋求最優解，是一種超啓發式算法，其中被稱爲個體的染色體（即研究對象的候選解）的集合能適應外部環境（研究對象的評價函數），並基於下面的規則在每代中生成新個體集合：
（1）越是適應性強的個體，其生存機率越高（自然選擇）；
（2）以原有個體爲基礎通過遺傳操作生成新個體（遺傳現象）。
其中（1）可看做是概率搜索法，（2）是經驗搜索法。

2.2 染色體設計（編碼和解碼）

揹包問題的候選解的編碼可採用一串長度爲物品個數的二進制列來表示，列中的每位表示物品的選擇情況，1表示物品被選擇，0表示物品未被選中。這樣一個一位二進制列即爲一個染色體（個體）。以下是算法的僞代碼，random[0,1]表示在[0,1]區間內隨機產生0或者1。

procedure:binary-string encoding
input:no.of item n
output:chromosome vk
begin
	for i=1 to n
		vk(i)<--random[0,1];
	output:chromosome vk;
end

編碼（種羣初始化）的python代碼如下：

#初始化種羣
def init(popsize,n): 
    population=[]
    for i in range(popsize):
        pop=''
        for j in range(n):
            pop=pop+str(np.random.randint(0,2))
        population.append(pop)    
    return population

將編碼表現爲候選解，就是使1對應的物體被選擇，因此解碼的僞代碼如下：

procedure:binary-string decoding
input:KP data set
		chromosome vk
output:profit f(vk),items in knapsack Sk
begin
	g<--0;
	f<--0;
	S<--0;
	for i=1 to n
		if vk=1 then
			if g+wi<=W then
				g<--g+wi;
				f<--f+ci;
				S<--S+{i};
			else
				break;
	end
	output:profit f(vk),items in knapsack Sk;
end

解碼和適應度函數計算的python代碼如下：

#解碼1
def decode1(x,n,w,c,W):
    s=[]#儲存被選擇物體的下標集合
    g=0
    f=0
    for i in range(n):
        if (x[i] == '1'):
            if g+w[i] <= W:
                g = g+w[i]
                f = f+c[i]
                s.append(i)
            else:
                break
    return f,s

#適應度函數1
def fitnessfun1(population,n,w,c,W):
    value=[]
    ss=[]
    for i in range(len(population)):
        [f,s]= decode1(population[i],n,w,c,W)
        value.append(f)
        ss.append(s)
    return value,ss

除了上述的解碼方式外，我額外提供一種帶懲罰函數的解碼以及適應度計算的方式，請讀者自行思考參閱。
這種帶有懲罰項的適應度計算方法如下：
$eval(v_k)=\left \{ \begin{matrix} \sum_{j\in S_k}c_j,\sum_{j\in S_k}w_j\leqslant W \\-M\sum_{j\in S_k}c_j,\sum_{j\in S_k}w_j> W \end{matrix} \right .$eval(vk​)={∑j∈Sk​​cj​,∑j∈Sk​​wj​⩽W−M∑j∈Sk​​cj​,∑j∈Sk​​wj​>W​
此方法的解碼和適應度函數計算的python代碼如下：

#解碼2
def decode2(x,n,w,c):
    s=[]#儲存被選擇物體的下標集合
    g=0
    f=0
    for i in range(n):
        if (x[i] == '1'):
            g = g+w[i]
            f = f+c[i]
            s.append(i)
    return g,f,s

#適應度函數2
def fitnessfun2(population,n,w,c,W,M):
    value=[]
    ss=[]
    for i in range(len(population)):
        [g,f,s]= decode2(population[i],n,w,c)
        if g>W:
            f = -M*f#懲罰
        value.append(f)
        ss.append(s)
    minvalue=min(value)
    value=[(i-minvalue+1) for i in value]
    return value,ss

2.3 遺傳操作

2.3.1輪盤賭模型

輪盤賭模型的基本原理是根據每個染色體的適值比例來確定該個體的選擇概率或生存概率。如下面的僞代碼那樣，個體適應度按比例轉換爲輪盤的面積並旋轉輪盤，最後選擇球落點位置所對應的個體。（此處輸入可以爲父代加子代，也可直接輸入子代，具體區別見後面的2.4算法流程）
這裏有一些參數說明：
l：染色體索引號
F：所有染色體的適值和
vk：第k個染色體
pk：第k個染色體的選擇概率
qk：前k個染色體的累計概率

procedure:roulette wheel selection
input:chromosome Pk,k=1,2,...,popsize+offsize
output:chromosome Pk in next generation
begin
	計算適應度函數值eval(Pk);
	計算累計適應度F=sum_{k=1->popsize+offsize}eval(Pk);
	計算選擇概率pk=eval(Pk)/F(k=1,2,...,popsize+offsize);
	計算累計概率qk=sum_{i=1->k}pi(k=1,2,...,popsize+offsize);
	for k=1 to popsize
		r<--random[0,1]
		if r<=q1 then
			Pl'<--Pk;
		else if qk-1<r<=qk then
			Pk'<--Pk;
	end
	output:chromosome Pk,k=1,2,...,popsize;
end

輪盤賭的python代碼如下;

#輪盤賭選擇
def roulettewheel(population,value,pop_num):
    fitness_sum=[]
    value_sum=sum(value)
    fitness=[i/value_sum for i in value]
    for i in range(len(population)):##
        if i==0:
            fitness_sum.append(fitness[i])
        else:
            fitness_sum.append(fitness_sum[i-1]+fitness[i])
    population_new=[]
    for j in range(pop_num):###
        r=np.random.uniform(0,1)
        for i in range(len(fitness_sum)):###
            if i==0:
                if r>=0 and r<=fitness_sum[i]:
                    population_new.append(population[i])
            else:
                if r>=fitness_sum[i-1] and r<=fitness_sum[i]:
                    population_new.append(population[i])
    return population_new

2.3.2兩點交叉

兩點交叉即在染色體中隨機產生兩個斷點，將這兩個斷點的中間部位進行交換的交叉方法。此算法父代的雙親個體有可能被多次選擇，也可能一次都未選到，即不滿足完備性（父代中的個體不是所有都被選擇）。此算法也可以讓父代的雙親個體都被選擇，即滿足完備性。以下的算法僞代碼爲前者，實現的python代碼爲後者。
參數說明：
pc：交叉概率
p：交叉點
L：爲染色體長度

input: pc, parent Pk, k=1, 2, ..., popsize 
output: offspring Ck
begin
	for k <-- 1 to popsize/2 do
		if pc≥random[0,1] then
			i<--0;
			j<--0;
			repeat
				i<--random[1,popsize];
				j<--random[1,popsize];
			until(i≠j)                
			p<--random[1,L-1]; 
			q<--random[p,L]; 
			Ck<--Pi[1:p]//Pj[p:q]//Pi[q:L];
			C(k+popsize/2)<--Pj[1:p]//Pi[p:q]//Pj[q:L];
		end
	end
	output:offspring Ck;
end

兩點交叉的python代碼如下：

#兩點交叉
def crossover(population_new,pc,ncross):
    a=int(len(population_new)/2)
    parents_one=population_new[:a]
    parents_two=population_new[a:]
    np.random.shuffle(parents_one)
    np.random.shuffle(parents_two)
    offspring=[]
    for i in range(a):
        r=np.random.uniform(0,1)
        if r<=pc:
            point1=np.random.randint(0,(len(parents_one[i])-1))
            point2=np.random.randint(point1,len(parents_one[i]))
            off_one=parents_one[i][:point1]+parents_two[i][point1:point2]+parents_one[i][point2:]
            off_two=parents_two[i][:point1]+parents_one[i][point1:point2]+parents_two[i][point2:]
            ncross = ncross+1
        else:
            off_one=parents_one[i]
            off_two=parents_two[i]
        offspring.append(off_one)
        offspring.append(off_two)
    return offspring

2.3.3反轉變異（單點變異）

實現單點變異算子的常用方法1爲每條染色體隨機生成一個數，此數指示該染色體是否需要變異。如果該染色體需要變異，爲其生成一個隨機變量，此隨機變量指示修改該染色體的哪一位。其算法的僞代碼如下：

procedure:mutation1
input: pm, parent Pk, k=1, 2, ..., popsize              // pm:變異概率
output: offspring Ck
begin
	for k <-- 1 to popsize do	                       
		if pm <-- random [0, 1] then
			point<--random[0,L-1]	          // L: 染色體長度
			if point = 0
				Pk <-- 1-Pk [ 0 ] // Pk[ 1: L ];	
			else
				Pk <-- Pk [1: point-1] // 1-Pk [ point ] // Pk[ point+1: L ];	
			end
		end
		Ck =Pk;
	end
   	output:offspring Ck;
end

反轉變異1的python代碼如下：

#單點變異1
def mutation1(offspring,pm,nmut):
    for i in range(len(offspring)):
        r=np.random.uniform(0,1)
        if r<=pm:
            point=np.random.randint(0,len(offspring[i]))
            if point==0:
                if offspring[i][point]=='1':
                    offspring[i]='0'+offspring[i][1:]
                else:
                    offspring[i]='1'+offspring[i][1:]
            else:
                if offspring[i][point]=='1':
                    offspring[i]=offspring[i][:(point-1)]+'0'+offspring[i][point:]
                else:
                    offspring[i]=offspring[i][:(point-1)]+'1'+offspring[i][point:]
            nmut = nmut+1
    return offspring

實現單點變異算子的常用方法2爲染色體數組中的每個位生成隨機變量。此隨機變量指示是否將修改特定位。其算法的僞代碼如下：

procedure:mutation2
input: pm, parent Pk, k=1, 2, ..., popsize              // pm:變異概率
output: offspring Ck
begin
	for k <-- 1 to popsize do	                       
		for j <-- 1 to L do(基因串長度）              // L: 染色體長度
			if pm <-- random [0, 1] then	          
				Pk <-- Pk [1: j-1] // 1-Pk [ j ] // Pk[ j+1: L ];	
			end
		end
		Ck =Pk;
	end
    output:offspring Ck;
end

反轉變異2的python代碼如下：

#單點變異2
def mutation2(offspring,pm,nmut):
    for i in range(len(offspring)):
        for j in range(len(offspring[i])):
            r=np.random.uniform(0,1)
            if r<=pm:
                if j==0:
                    if offspring[i][j]=='1':
                        offspring[i]='0'+offspring[i][1:]
                    else:
                        offspring[i]='1'+offspring[i][1:]
                else:
                    if offspring[i][j]=='1':
                        offspring[i]=offspring[i][:(j-1)]+'0'+offspring[i][j:]
                    else:
                        offspring[i]=offspring[i][:(j-1)]+'1'+offspring[i][j:]
                nmut = nmut+1
    return offspring

2.4 算法流程

應用經典遺傳算法對01揹包問題求解的總體流程僞代碼如下：（此處由於對輪盤賭操作的輸入種羣不同，因此列出兩種python代碼，請讀者加以區別。）
GA的算法流程爲僞代碼如下：

procedure:GA for Knapsack Problem(KP)
input: KP data set, GA parameters
output: the best solution
begin
	t<-- 0	                       // t:遺傳代數
	initialize P(t) by encoding routine;         //P(t):染色體種羣 
	fitness eval(P) by decoding routine;
	while (not termination condition) do
		crossover P(t) to yield C(t);     //C(t):offspring
		mutation   P(t) to yield C(t);
		fitness eval(C) by decoding routine;
		select P(t+1) from P(t) and C(t);
		t<--t+1;
	end
    output:the best solution;
end

GA的python代碼1如下（輪盤賭的輸入爲父代加子代）：

import numpy as np
import matplotlib.pyplot as plt
#主程序----------------------------------------------------------------------------------------------------------------------------------
#參數設置-----------------------------------------------------------------------
gen=1000#迭代次數
pc=0.25#交叉概率
pm=0.02#變異概率
popsize=10#種羣大小
n = 20#物品數,即染色體長度n
w = [2,5,18,3,2,5,10,4,8,12,5,10,7,15,11,2,8,10,5,9]#每個物品的重量列表
c = [5,10,12,4,3,11,13,10,7,15,8,19,1,17,12,9,15,20,2,6]#每個物品的代價列表
W = 40#揹包容量
M = 5#懲罰值
fun = 1#1-第一種解碼方式，2-第二種解碼方式（懲罰項） 
#初始化-------------------------------------------------------------------------
#初始化種羣（編碼）
population=init(popsize,n)
#適應度評價（解碼）
if fun==1:
    value,s = fitnessfun1(population,n,w,c,W)
else:
    value,s = fitnessfun2(population,n,w,c,W,M)
#初始化交叉個數
ncross=0
#初始化變異個數
nmut=0
#儲存每代種羣的最優值及其對應的個體
t=[]
best_ind=[]
last=[]#儲存最後一代個體的適應度值
realvalue=[]#儲存最後一代解碼後的值
#循環---------------------------------------------------------------------------
for i in range(gen):
    print("迭代次數：")
    print(i)
    #交叉
    offspring_c=crossover(population,pc,ncross)
    #變異
    #offspring_m=mutation1(offspring,pm,nmut)
    offspring_m=mutation2(offspring_c,pm,nmut)
    mixpopulation=population+offspring_m
    #適應度函數計算
    if fun==1:
        value,s = fitnessfun1(mixpopulation,n,w,c,W)
    else:
        value,s = fitnessfun2(mixpopulation,n,w,c,W,M)
    #輪盤賭選擇
    population=roulettewheel(mixpopulation,value,popsize)
    #儲存當代的最優解
    result=[]
    if i==gen-1:
        if fun==1:
            value1,s1 = fitnessfun1(population,n,w,c,W)
            realvalue=s1
            result=value1
            last=value1
        else:
            for j in range(len(population)):
                g1,f1,s1 = decode2(population[j],n,w,c)
                result.append(f1)
                realvalue.append(s1)
            last=result
    else:
        if fun==1:
            value1,s1 = fitnessfun1(population,n,w,c,W)
            result=value1
        else:
            for j in range(len(population)):
                g1,f1,s1 = decode2(population[j],n,w,c)
                result.append(f1)
    maxre=max(result)
    h=result.index(max(result))
    #將每代的最優解加入結果種羣
    t.append(maxre)
    best_ind.append(population[h])

#輸出結果-----------------------------------------------------------------------
if fun==1:
    best_value=max(t)
    hh=t.index(max(t))
    f2,s2 = decode1(best_ind[hh],n,w,c,W)
    print("最優組合爲：")
    print(s2)
    print("最優解爲：")
    print(f2)
    print("最優解出現的代數：")
    print(hh)
    #畫出收斂曲線
    plt.plot(t)
    plt.title('The curve of the optimal function value of each generation with the number of iterations',color='#123456')
    plt.xlabel('the number of iterations')
    plt.ylabel('the optimal function value of each generation')
else:
    best_value=max(result)
    hh=result.index(max(result))
    s2 = realvalue[hh]
    print("最優組合爲：")
    print(s2)
    print("最優解爲：")
    print(f2)

GA的python代碼2如下（輪盤賭的輸入爲子代）：

import numpy as np
import matplotlib.pyplot as plt

#主程序----------------------------------------------------------------------------------------------------------------------------------
#參數設置-----------------------------------------------------------------------
gen=1000#迭代次數
pc=0.25#交叉概率
pm=0.02#變異概率
popsize=20#種羣大小
n = 20#物品數,即染色體長度n
w = [2,5,18,3,2,5,10,4,8,12,5,10,7,15,11,2,8,10,5,9]#每個物品的重量列表
c = [5,10,12,4,3,11,13,10,7,15,8,19,1,17,12,9,15,20,2,6]#每個物品的代價列表
W = 40#揹包容量
M = 5#懲罰值
fun = 1#1-第一種解碼方式，2-第二種解碼方式（懲罰項） 
#初始化-------------------------------------------------------------------------
#初始化種羣（編碼）
population=init(popsize,n)
#適應度評價（解碼）
if fun==1:
    value,s = fitnessfun1(population,n,w,c,W)
else:
    value,s = fitnessfun2(population,n,w,c,W,M)
#初始化交叉個數
ncross=0
#初始化變異個數
nmut=0
#儲存每代種羣的最優值及其對應的個體
t=[]
best_ind=[]
last=[]#儲存最後一代個體的適應度值
realvalue=[]#儲存最後一代解碼後的值
#循環---------------------------------------------------------------------------
for i in range(gen):
    print("迭代次數：")
    print(i)
    #適應度函數計算
    if fun==1:
        value,s = fitnessfun1(population,n,w,c,W)
    else:
        value,s = fitnessfun2(population,n,w,c,W,M)
    #輪盤賭選擇
    population=roulettewheel(population,value,popsize)
    #交叉
    offspring_c=crossover(population,pc,ncross)
    #變異
    #offspring_m=mutation1(offspring,pm,nmut)
    offspring_m=mutation2(offspring_c,pm,nmut)
    population=offspring_m

    #儲存當代的最優解
    result=[]
    if i==gen-1:
        if fun==1:
            value1,s1 = fitnessfun1(population,n,w,c,W)
            realvalue=s1
            result=value1
            last=value1
        else:
            for j in range(len(population)):
                g1,f1,s1 = decode2(population[j],n,w,c)
                result.append(f1)
                realvalue.append(s1)
            last=result
    else:
        if fun==1:
            value1,s1 = fitnessfun1(population,n,w,c,W)
            result=value1
        else:
            for j in range(len(population)):
                g1,f1,s1 = decode2(population[j],n,w,c)
                result.append(f1)
    maxre=max(result)
    h=result.index(max(result))
    #將每代的最優解加入結果種羣
    t.append(maxre)
    best_ind.append(population[h])

#輸出結果-----------------------------------------------------------------------
if fun==1:
    best_value=max(t)
    hh=t.index(max(t))
    f2,s2 = decode1(best_ind[hh],n,w,c,W)
    print("最有組合爲：")
    print(s2)
    print("最優解爲：")
    print(f2)
    print("最優解出現的代數：")
    print(hh)
    #畫出收斂曲線
    plt.plot(t)
    plt.title('The curve of the optimal function value of each generation with the number of iterations',color='#123456')
    plt.xlabel('the number of iterations')
    plt.ylabel('the optimal function value of each generation')
else:
    best_value=max(result)
    hh=result.index(max(result))
    s2 = realvalue[hh]
    print("最優組合爲：")
    print(s2)
    print("最優解爲：")
    print(f2)

3. 實驗結果

本實驗參數設置如下：
1.gen=1000#迭代次數
2.pc=0.25#交叉概率
3.pm=0.02#變異概率
4.popsize=20#種羣大小
5.n = 20#物品數,即染色體長度n，要與6，7的列表維度保持一致
6.w = [2,5,18,3,2,5,10,4,8,12,5,10,7,15,11,2,8,10,5,9]#每個物品的重量列表
7.c = [5,10,12,4,3,11,13,10,7,15,8,19,1,17,12,9,15,20,2,6]#每個物品的代價列表
8.W = 40#揹包容量
9.M = 5#懲罰值
根據上述參數設置和原理分析，我們編程得到每代最優的函數值隨進化代數的變化如下圖所示。由於兩種不同的解碼和適應度計算的方式，最後二者的結果呈現方式不同。換句話說，這個圖是第一種解碼和適應度計算的方式得到的圖，請思考爲什麼第二種方法沒有畫這個曲線？（提示：帶懲罰值的方法每代最優解的適應度值並非揹包問題的函數值。）

完整的python代碼如下：https://download.csdn.net/download/qq_40434430/10743582 由於作者水平有限，以上內容難免有錯誤之處，如果讀者有所發現，請在下面留言告知我，鄙人將不勝感激！ [1]:(日)玄光男, 林林等著. 網絡模型與多目標遺傳算法[M]. 於歆傑, 梁承姬譯. 北京: 清華大學出版社, 2017