sopt:一個簡單的python最優化庫

引言

    最近有些朋友總來問我有關遺傳算法的東西，我是在大學搞數學建模的時候接觸過一些最優化和進化算法方面的東西，以前也寫過幾篇博客記錄過,比如[遺傳算法的C語言實現(一):以非線性函數求極值為例](https://www.cnblogs.com/lyrichu/p/6152897.html)和[C語言實現粒子群算法(PSO)一](https://www.cnblogs.com/lyrichu/p/6151272.html)等，如果對原理有興趣的話可以去我的博客具體查看:[Lyrichu's Blog](https://www.cnblogs.com/lyrichu)。所以突發奇想，干脆把以前寫的一些進化算法比如遺傳算法(GA),粒子群算法(PSO),模擬退火算法(SA)以及最近看的基於梯度的一些優化算法比如Gradient Descent,SGD,Momentum等整理一下，寫成一個python庫，方便那些有需要的朋友使用。斷斷續續花了幾天的時間，初步完成了一些基本功能，這里簡單介紹一下。

sopt 簡介

    **sopt**是simple optimization的簡稱，目前我已經將代碼托管到pypi,地址是[sopt](https://pypi.org/project/sopt/),可以直接通過pip命令下載安裝使用,由於項目只依賴numpy,所以在windows和linux環境下安裝都很方便，直接`pip install sopt`即可。項目的github地址是[sopt](https://github.com/Lyrichu/sopt)。目前sopt包含的優化方法如下: - 遺傳算法(Genetic Algorithm,GA) - 粒子群算法(Particle Swarm Optimization,PSO) - 模擬退火算法(Simulated Anealing,SA) - 隨機游走算法(Random Walk): - 梯度下降法(Gradient Descent,GD) - 動量優化算法(Momentum) - 自適應梯度算法(AdaGrad) - RMSProp - Adam 由於只是一個初步的版本，后續如果有時間的話，會加上更多的優化算法進去。目前所有的優化算法暫時只支持**連續實函數**的優化;除了基於梯度的幾個優化算法,GA、PSO、SA以及Random Walk都同時支持**無約束優化**，**線性約束優化**以及**非線性約束優化**。具體我會在下面詳細說明。

sopt 使用詳解以及實例演示

1.SGA 使用

    SGA是Simple Genetic Algorithm的簡稱，是最基本的遺傳算法。其編碼方式采用**二進制編碼**,選擇方法采用**輪盤賭法**,交叉方法采用**單點交叉**,變異方式采用**均勻變異**,默認是求函數的**最小值**。下面是sopt中**SGA**的一個簡單使用實例: ```python from sopt.SGA import SGA from math import sin def func1(x): return (x[0]-1)**2 + (sin(x[1])-0.5)**4 + 2 if __name__ == '__main__': sga = SGA.SGA(func = func1,func_type = 'min',variables_num = 2, lower_bound = 0,upper_bound = 2,generations = 20, binary_code_length = 10) # run SGA sga.run() # show the SGA optimization result in figure sga.save_plot() # print the result sga.show_result() ``` 運行結果如下: ``` -------------------- SGA config is: -------------------- lower_bound:[0, 0] generations:20 cross_rate:0.7 variables_num:2 mutation_rate:0.1 func_type:min upper_bound:[2, 2] population_size:100 func: binary_code_length:10 -------------------- SGA caculation result is: -------------------- global best generation index/total generations:3/20 global best point:[1.00488759 0.45356794] global best target:2.00003849823336 ``` 用圖像展示為圖1所示: 圖1 SGA 運行結果 上面定義的目標函數為$f(x_1,x_2)=(x_1-1)^2+(sin(x_2)-0.5)^4+2$,其中$0 2. GA 使用     GA是相比SGA更加一般的遺傳算法實現，其對於編碼方式，選擇方式，交叉方式以變異方式等會有更多的支持。首先還是看一個例子: ```python from sopt.GA.GA import GA from sopt.util.functions import * from sopt.util.ga_config import * from sopt.util.constraints import * class TestGA: def __init__(self): self.func = quadratic11 self.func_type = quadratic11_func_type self.variables_num = quadratic11_variables_num self.lower_bound = quadratic11_lower_bound self.upper_bound = quadratic11_upper_bound self.cross_rate = 0.8 self.mutation_rate = 0.05 self.generations = 200 self.population_size = 100 self.binary_code_length = 20 self.cross_rate_exp = 1 self.mutation_rate_exp = 1 self.code_type = code_type.binary self.cross_code = False self.select_method = select_method.proportion self.rank_select_probs = None self.tournament_num = 2 self.cross_method = cross_method.uniform self.arithmetic_cross_alpha = 0.1 self.arithmetic_cross_exp = 1 self.mutation_method = mutation_method.uniform self.none_uniform_mutation_rate = 1 #self.complex_constraints = [constraints1,constraints2,constraints3] self.complex_constraints = None self.complex_constraints_method = complex_constraints_method.penalty self.complex_constraints_C = 1e6 self.M = 1e8 self.GA = GA(**self.__dict__) def test(self): start_time = time() self.GA.run() print("GA costs %.4f seconds!" % (time()-start_time)) self.GA.save_plot() self.GA.show_result() if __name__ == '__main__': TestGA().test() ``` 上面代碼的運行結果為: ``` GA costs 6.8320 seconds! -------------------- GA config is: -------------------- func: code_type:binary complex_constraints:None global_generations_step:200 cross_method:uniform mutation_method:uniform cross_rate:0.8 lower_bound:[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] tournament_num:2 variables_num:11 complex_constraints_method:penalty none_uniform_mutation_rate:1 population_size:100 mutation_rate:0.05 generations:200 arithmetic_cross_alpha:0.1 func_type:min mutation_rate_exp:1 cross_rate_exp:1 arithmetic_cross_exp:1 M:100000000.0 select_method:proportion upper_bound:[11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11] cross_code:False binary_code_length:20 complex_constraints_C:1000000.0 -------------------- GA caculation result is: -------------------- global best target generation index/total generations:149/200 global best point:[ 1.07431037 2.41401426 2.4807906 4.36291634 4.90653029 6.13753427 6.58147963 7.7370479 9.42957347 10.46616122 10.87151134] global best target:2.2685208668743204 ```     其中`sopt.util.functions`預定義了一些測試函數,`quadratic11`是一個11元變量的二次函數，函數原型為:$quadratic11(x_1,...,x_{11})=(x_1-1)^2 +(x_2-2)^2 + ... +(x_{11}-11)^2 + 1$,其中$1 \le x_1,...,x_{11} \le 11$,函數的最小值點在(1,2,...,11)處取得，最小值為1。另外還定義了其他幾個測試函數為: - quadratic50:和quadratic11類似,只是變量個數變為50,取值范圍變為1-50; - quadratic100:和quadratic11類似,只是變量個數變為100,取值范圍變為1-100; - quadratic500:和quadratic11類似,只是變量個數變為500,取值范圍變為1-500; - quadratic1000:和quadratic1000類似,只是變量個數變為1000,取值范圍變為1-1000; - Rosenbrock:函數原型為:$Rosenbrock(x_1,x_2)=100(x_1^2-x_2)^2+(1-x_1)^2$,其中$-2.048 \le x_1,x_2 \le 2.048$,這個函數有很多局部極小值點，函數在(-2.048,-2.048)處取得最大值; - shubert2函數:函數定義為,$shubert2(x_1,x_2)=\prod_{i=1}^{2}(\sum_{j=1}^{5}(jcos((j+1)*x_i+j))$,其中$-10 \le x_1,x_2 \le 10$,函數有760個局部極小值點,全局極小值為-186.731,這里為了將適應度函數變換為正，我們在原函數的基礎上加上1000,這樣函數的全局極小值點就變為813.269; - shubert函數，將shubert2函數中的兩個變量拓展為n個變量($\prod_{i=1}^{n}$)就可以得到一般的shubert函數了。     對於每一個優化方法，我們都預定了一個形如`xx_config`的模塊，其中定義了該優化方法的一些常用默認參數，比如`ga_config`中就定義了一些GA的一些常用優化參數，`ga_config.basic_config`定義了一些基礎參數設置，比如`basic_config.generations`是一個默認進化代數,`basic_config.mutation_rate`是默認的變異參數；而`ga_config.cross_method`則預定義了所有支持的交叉方法,比如`cross_method.uniform`表示均勻交叉,`cross_method.one_point`表示單點交叉等;`ga_config.mutation_method`等也是類似的。有了這些預定義變量，可以免去我們手動輸入很多參數取值以及傳入方法字符串的麻煩(有時候可能會寫錯)。     觀察上面的運行結果，我們發現對於quadratic11函數，最終找到的全局極小值為2.2685208668743204,和真實的全局極小值點1已經比較接近了；運行耗時6秒多，似乎有些長，這是因為我們種群數目設為100，進化代數設為200,比較大，而且又是采用二進制20位編碼，再加上python腳本語言的運行效率問題，所以稍慢。為了做一個對比,這里我們將目標函數從quadratic11變為Rosenbrock函數，其他參數設置保持不變，得到的結果如下: ``` GA costs 1.7245 seconds! -------------------- GA config is: -------------------- population_size:100 lower_bound:[-2.048, -2.048] mutation_rate_exp:1 select_method:proportion code_type:binary global_generations_step:200 generations:200 mutation_method:uniform complex_constraints_method:penalty binary_code_length:20 cross_method:uniform arithmetic_cross_alpha:0.1 func: upper_bound:[2.048, 2.048] cross_code:False mutation_rate:0.05 cross_rate_exp:1 complex_constraints_C:1000000.0 cross_rate:0.8 variables_num:2 M:100000000.0 complex_constraints:None none_uniform_mutation_rate:1 tournament_num:2 arithmetic_cross_exp:1 func_type:max -------------------- GA caculation result is: -------------------- global best target generation index/total generations:75/200 global best point:[-2.04776953 -2.04537109] global best target:3901.4655271502425 ``` 圖2是每代最優值的計算結果: 圖2 GA Rosenbrock 函數運行200代尋優結果 替換成Rosenbrock函數之后，函數的運行時間從6秒多減少到1秒多(函數變量個數明顯減少了),這說明GA的運行時間與目標函數的變量個數是顯著相關的。最后找到的全局極大值點為(-2.04776953,-2.04537109)和真實全局極大值點(-2.048,-2.048)已經很接近了。     SGA類是GA類的父類，所以GA類和SGA類有一些公共的屬性，比如`generations`,`population_size`,`func_type`等。和SGA一樣的參數這里就不再列舉了，如下是GA特有的一些參數: - cross_rate_exp:cross_rate指數遞增的參數取值,即變異率按照$r_t=r_0 \beta^t$,這里$r_t$表示t次迭代的交叉率,$r_0$是初始的交叉率，$\beta$就是這里的`cross_rate_exp`,默認取值為1,一般設置為一個比1稍大的數字，比如1.0001等; - mutation_rate_exp:mutation_rate指數遞增參數取值,具體含義和`cross_rate_exp`類似; - code_type:采用什么樣的編碼方式,有三種取值:'binary'表示二進制編碼(默認),'gray'表示采用格雷編碼,'real'表示采用實數編碼。其中格雷編碼是一種改進的二進制編碼，其優點在於進行交叉、變異等操作時，改變了染色體基因的某幾個位置，二進制編碼對應的實數值可能會發生非常大的變化，而格雷編碼一般不會。二進制編碼與格雷編碼的轉換可以參考附錄; - cross_code:這是一個布爾值,表示是否對二進制編碼或格雷碼采用交叉編碼的方式。具體說明參考附錄; - select_method:選擇操作方法，有如下的6種取值: 1. 'proportion':表示比例選擇(輪盤賭法) 2. 'keep_best':表示保留最優個體 3. 'determinate_sampling':表示確定性采樣 4. 'rssr':remainder stochastic sampling with replacement的縮寫，表示無放回余數隨機選擇 5. 'rank':表示排序選擇 6. 'stochastic_tournament':表示隨機錦標賽法(以上6種選擇方式詳細說明請參考附錄) - rank_select_probs:僅當select_method取值為'rank'時，該參數才起作用，表示采用排序選擇時，適應度從低到高每個個體被選擇的概率，默認為None,采用$p_i=\frac{i}{\sum_{j=1}^{n}j},i=1,2,...,n$的方式計算，其中$n$表示種群個數,$p_i$表示第$i$個個體被選擇的概率，所有概率之和為1，如果取值不為None,你應該傳入一個size為`population_size`的ndarray,所有元素按照遞增排序，數組和為1; - tournament_num:如果select_method取值為'stochastic_tournament',該值表示錦標賽競爭者數量; - cross_method:交叉方法,有如下4種取值: 1. 'one_point'：表示單點交叉 2. 'two_point':表示雙點交叉 3. 'uniform':表示均勻交叉 4. 'arithmetic':表示算術交叉,只有當編碼采用實數編碼時，才能使用算術交叉 - arithmetic_cross_alpha:算術交叉的系數，算術交叉的公式為:$x_{new1}=\alpha x_1 + (1-\alpha)x_2,x_{new2}=\alpha x_2 + (1-\alpha)x_1$,這里的$\alpha$即表示`arithmetic_cross_alpha`,其默認值為0.1; - 'arithmetic_cross_exp':算術交叉系數按照指數遞增，即$\alpha_t = \alpha r^t$,這里$\alpha_t$表示第t代的算術交叉系數,$\alpha$是初始交叉系數,$r$即是這里的`arithmetic_cross_exp`,默認取值為1，一般取一個比1稍大的數，比如1.0001,t是進化代數; - 'mutation_method':變異方法，有如下5種取值: 1. simple:即簡單變異，隨機選擇一個染色體上的基因，按照變異概率進行變異 2. uniform:即均勻變異,每個染色體上的每個基因都按照變異概率進行變異 3. boundary:邊界變異,每個基因進行變異以后的取值只能去邊界值,比如說各以0.5的概率取得其上界和下界，這種變異方式僅適用於實數編碼的情況，一般在目標函數的最優點靠近邊界時使用 4. none_uniform:非均勻變異，我們不是取均勻分布的隨機數去替換原來的基因，而是在原來基因的基礎上做一點微小的隨機擾動，擾動以后的結果作為新的基因值，具體來說，我們采用的是這樣的變異方式:1)if random(0,1) = 0,$x_k^{'} = x_k + \bigtriangleup (t,U_{max}^{k}-x_k)$;2)if random(0,1) = 1,$x_k{'} = x_k - \bigtriangleup (t,x_k -U_{min}{k})$ 其中 $\bigtriangleup(t,y)=y(1-r^{(1-t/T)b})$,$r$是0-1均勻分布的一個隨機數，$T$是最大進化代數，$b$是一個系統參數，表示隨機擾動對於進化代數的依賴長度，默認值為1 5. gaussian:高斯變異,即使用高斯分布去替代均勻分布來進行變異,由高斯分布的特點可以知道，這種變異方式也是在原個體區域附近的某個局部區域進行重點搜索，這里高斯分布的均值$\mu$定義為$\mu = \frac{U_{min}^{k}+U_{max}^{k}}{2}$,標准差$\sigma=\frac{U_{max}^{k}-U_{min}^{k}}{6}$ - none_uniform_mutation_rate:即均勻分布`none_uniform`中定義的系統參數$b$; - complex_constraints:復雜約束,默認值為None,即沒有復雜約束,只有簡單的邊界約束,如果有復雜約束，其取值應該是$\[func_1,func_2,...,func_n\]$, 其中$func_i$表示第$i$個復雜約束函數名，比如對於一個復雜約束函數$func_1:x_1^2+x_2^2 < 3$,其復雜約束函數應該這樣定義: ``` def func1(x): x1 = x[0] x2 = x[1] return x1**2 + x2**2 - 3 ``` - complex_constraints_method:復雜約束求解的方法,默認是`penalty`即懲罰函數法，暫時不支持其他的求解方式; - complex_constraints_C:采用`penalty`求解復雜約束的系數$C$,比如對於某一個約束$x_1^2+x_2^2 < 3$,GA在求解過程中，違反了該約束,即解滿足$x_1^2+x_2^2 \ge 3$,那么我們對目標函數增加一個懲罰項: $C(x_1^2+x_2^2-3)$,$C$一般取一個很大的正數，默認值為$10^6$。     GA類的參數很多可以調節，所以相比之下使用起來更加麻煩，特別是對於復雜函數的尋優，默認參數未必是最好的，可能需要你根據目標函數的特點自行嘗試，選擇最優的參數組合。

3. PSO 使用

    PSO算法的原理可以參考我之前的博客[C語言實現粒子群算法(PSO)一](http://www.cnblogs.com/lyrichu/p/6151272.html)和[C語言實現粒子群算法（PSO）二](http://www.cnblogs.com/lyrichu/p/6151293.html),這里不再贅述。還是先看一個實例代碼: ```python from time import time from sopt.util.functions import * from sopt.util.pso_config import * from sopt.PSO.PSO import PSO from sopt.util.constraints import * class TestPSO: def __init__(self): self.func = quadratic11 self.func_type = quadratic11_func_type self.variables_num = quadratic11_variables_num self.lower_bound = quadratic11_lower_bound self.upper_bound = quadratic11_upper_bound self.c1 = basic_config.c1 self.c2 = basic_config.c2 self.generations = 200 self.population_size = 100 self.vmax = 1 self.vmin = -1 self.w = 1 self.w_start = 0.9 self.w_end = 0.4 self.w_method = pso_w_method.linear_decrease #self.complex_constraints = [constraints1,constraints2,constraints3] self.complex_constraints = None self.complex_constraints_method = complex_constraints_method.loop self.PSO = PSO(**self.__dict__) def test(self): start_time = time() self.PSO.run() print("PSO costs %.4f seconds!" %(time()-start_time)) self.PSO.save_plot() self.PSO.show_result() if __name__ == '__main__': TestPSO().test() ``` 運行結果為: ``` PSO costs 1.1731 seconds! -------------------- PSO config is: -------------------- complex_constraints_method:loop c1:1.49445 lower_bound:[1 1 1 1 1 1 1 1 1 1 1] w_end:0.4 w_method:linear_decrease complex_constraints:None func: upper_bound:[11 11 11 11 11 11 11 11 11 11 11] generations:200 func_type:min w:1 c2:1.49445 w_start:0.9 population_size:100 vmin:-1 vmax:1 variables_num:11 -------------------- PSO calculation result is: -------------------- global best generation index/total generations: 198/200 global best point: [ 1. 1.99999999 2.99999999 4. 5. 6. 7.00000001 7.99999999 9.00000001 10.00000001 11. ] global best target: 1.0 ``` 圖3 PSO 求解 quadratic11 200代運行結果 上面的代碼意圖應該是非常明顯的,目標函數是$quadratic11$,最終求得的最小值點幾乎就是全局最小值點。下面是PSO類中所有參數的具體定義: - variables_num:變量個數(必填); - lower_bound:一個int或者float的數字，或者是長度為`variables_num`的ndarray,表示每個變量的下界,如果是一個數字的話，我們認為所有的下界都是一樣的(必填); - upper_bound:一個int或者float的數字,或者是長度為`variables_num`的ndarray,表示每個變量的上界,如果是一個數字的話，我們認為所有的上界都是一樣的(必填); - func:目標函數名(必填); - func_type:函數優化類型,求最小值取值為'min',求最大值取值為'max',默認是'min'; - c1:PSO 參數c1,默認值為1.49445; - c2:PSO 參數c2,默認值為1.49445; - generations:進化代數,默認值為100; - population_size:種群數量,默認為50; - vmax:粒子最大速度,默認值為1; - vmin:粒子最小速度,默認值為-1; - w:粒子慣性權重,默認為1; - w_start:如果不是使用恆定的慣性權重話,比如使用權重遞減策略,w_start表示初始權重; - w_end:相應的,w_end表示末尾權重; - w_method:權重遞減的方式,有如下5種取值: 1. 'constant':權重恆定 2. 'linear_decrease':線性遞減,即$w_t = w_{end} + (w_{start} - w_{end})\frac{(T-t)}{T}$,其中$T$表示最大進化代數,$w_t$表示第$t$代權重 3. 'square1_decrease':第一種平方遞減,即:$w_t = w_{start} - (w_{start}-w_{end})(\frac{t}{T})^2$ 4. 'square2_decrease':第二種平方遞減,即：$w_t = w_{start} - (w_{start}-w_{end})(\frac{2t}{T}-(\frac{t}{T})^2)$ 5. 'exp_decrease':指數遞減,即:$w_t = w_{end}(\frac{w_{start}}{w_{end}})^{\frac{1}{1+\frac{10t}{T}}}$ - complex_constraints:復雜約束，默認為None,具體含義參考GA類的complex_constraints - complex_constraints_method:求解復雜約束的方法，默認為'loop',即如果解不滿足復雜約束，則再次隨機產生解，直到滿足約束,暫時不支持其他的求解方式。

4. SA 使用

    SA,即模擬退火算法，是基於概率的一種尋優方法，具體原理可以參考我的博客[模擬退火算法（SA）求解TSP 問題（C語言實現）](http://www.cnblogs.com/lyrichu/p/6688459.html)。sopt中SA的使用實例如下: ```python from time import time from sopt.util.functions import * from sopt.Optimizers.SA import SA from sopt.util.sa_config import * from sopt.util.constraints import * class TestSA: def __init__(self): self.func = Rosenbrock self.func_type = Rosenbrock_func_type self.variables_num = Rosenbrock_variables_num self.lower_bound = Rosenbrock_lower_bound self.upper_bound = Rosenbrock_upper_bound self.T_start = 100 self.T_end = 1e-6 self.q = 0.9 self.L = 100 self.init_pos = None #self.complex_constraints = [constraints1,constraints2,constraints3] self.complex_constraints_method = complex_constraints_method.loop self.SA = SA(**self.__dict__) def test(self): start_time = time() self.SA.run() print("SA costs %.4f seconds!" %(time()-start_time)) self.SA.save_plot() self.SA.show_result() if __name__ == '__main__': TestSA().test() ``` 運行結果如下: ``` SA costs 0.2039 seconds! -------------------- SA config is: -------------------- func_type:max complex_constraints:None q:0.9 complex_constraints_method:loop T_start:100 steps:17500 T_end:1e-06 L:100 func: variables_num:2 lower_bound:[-2.048 -2.048] init_pos:[-2.03887265 -2.02503927] upper_bound:[2.048 2.048] -------------------- SA calculation result is: -------------------- global best generation index/total generations:2126/17500 global best point: [-2.03887265 -2.02503927] global best target: 3830.997799328349 ``` 圖4 SA求解Rosenbrock函數 SA類的具體參數含義如下: - variables_num:變量個數(必填); - lower_bound:一個int或者float的數字，或者是長度為`variables_num`的ndarray,表示每個變量的下界,如果是一個數字的話，我們認為所有的下界都是一樣的(必填); - upper_bound:一個int或者float的數字,或者是長度為`variables_num`的ndarray,表示每個變量的上界,如果是一個數字的話，我們認為所有的上界都是一樣的(必填); - func:目標函數名(必填); - func_type:函數優化類型,求最小值取值為'min',求最大值取值為'max',默認是'min'; - T_start:初始溫度,默認值為100; - T_end:末尾溫度,默認值為$10^{-6}$; - q:退火系數,默認值為0.9,一般取一個在[0.9,1)之間的數字; - L:每個溫度時的迭代次數，即鏈長,默認值為100; - init_pos:初始解的取值，默認值為None,此時初始解是隨機生成的,或者你也可以指定一個用ndarray表示的初始解; - complex_constraints:復雜約束，默認為None,表示沒有約束，具體定義同GA 的 complex_constraints; - complex_constraints_method:求解復雜約束的方法，默認為'loop',即如果解不滿足復雜約束，則再次隨機產生解，直到滿足約束,暫時不支持其他的求解方式。

5. Random Walk 使用

    Random Walk 即隨機游走算法，是一種全局隨機搜索算法，具體原理可以參考我的博客[介紹一個全局最優化的方法：隨機游走算法(Random Walk)](http://www.cnblogs.com/lyrichu/p/7209529.html)。sopt中Random Walk的使用實例如下: ```python from time import time from sopt.util.functions import * from sopt.util.constraints import * from sopt.util.random_walk_config import * from sopt.Optimizers.RandomWalk import RandomWalk class TestRandomWalk: def __init__(self): self.func = quadratic50 self.func_type = quadratic50_func_type self.variables_num = quadratic50_variables_num self.lower_bound = quadratic50_lower_bound self.upper_bound = quadratic50_upper_bound self.generations = 100 self.init_step = 10 self.eps = 1e-4 self.vectors_num = 10 self.init_pos = None # self.complex_constraints = [constraints1,constraints2,constraints3] self.complex_constraints = None self.complex_constraints_method = complex_constraints_method.loop self.RandomWalk = RandomWalk(**self.__dict__) def test(self): start_time = time() self.RandomWalk.random_walk() print("random walk costs %.4f seconds!" %(time() - start_time)) self.RandomWalk.save_plot() self.RandomWalk.show_result() if __name__ == '__main__': TestRandomWalk().test() ``` 運行結果為: ``` Finish 1 random walk! Finish 2 random walk! Finish 3 random walk! Finish 4 random walk! Finish 5 random walk! Finish 6 random walk! Finish 7 random walk! Finish 8 random walk! Finish 9 random walk! Finish 10 random walk! Finish 11 random walk! Finish 12 random walk! Finish 13 random walk! Finish 14 random walk! Finish 15 random walk! Finish 16 random walk! Finish 17 random walk! random walk costs 1.0647 seconds! -------------------- random walk config is: -------------------- init_step:10 eps:0.0001 generations_nums:9042 lower_bound:[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1] complex_constraints_method:loop walk_nums:17 complex_constraints:None vectors_num:10 upper_bound:[50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50] variables_num:50 func: generations:100 func_type:min -------------------- random walk caculation result is: -------------------- global best generation index/total generations:8942/9042 global best point is: [ 1.00004803 1.99999419 3.00000569 3.99998558 5.00002455 5.99999255 6.99992476 7.99992864 9.00000401 9.99994717 10.99998155 12.00002429 13.0000035 13.99998567 15.00000421 16.00001454 16.99997252 17.99998041 19.00002491 20.00003141 21.00004182 21.99998565 22.99997668 23.99999821 24.99995881 25.99999359 27.00000443 28.00005117 28.99998132 30.00004136 31.00002021 32.00000616 33.00000678 34.00005423 35.00001799 36.00000051 37.00002749 38.00000203 39.00007087 39.9999964 41.00004432 42.0000158 42.99992991 43.99995352 44.99997267 46.00003533 46.9999834 47.99996778 49.00002904 50. ] global best target is: 1.0000000528527013 ``` 圖5 Random Walk 求解quadratic50 經過實驗發現,Random Walk 具有非常強的全局尋優能力，對於quadratic50這種具有50個變量的復雜目標函數，它也可以很快找到其全局最優點，而且運行速度也很快。RandomWalk類的具體參數含義如下: - variables_num:變量個數(必填); - lower_bound:一個int或者float的數字，或者是長度為`variables_num`的ndarray,表示每個變量的下界,如果是一個數字的話，我們認為所有的下界都是一樣的(必填); - upper_bound:一個int或者float的數字,或者是長度為`variables_num`的ndarray,表示每個變量的上界,如果是一個數字的話，我們認為所有的上界都是一樣的(必填); - func:目標函數名(必填); - func_type:函數優化類型,求最小值取值為'min',求最大值取值為'max',默認是'min'; - generations:每個step的最大迭代次數,默認是100; - init_step:初始步長(step)，默認值為10.0; - eps:終止迭代的步長,默認為$10^{-4}$; - vectors_num:使用 improved_random_walk時隨機產生向量的個數,默認為10; - init_pos:初始解的取值，默認值為None,此時初始解是隨機生成的,或者你也可以指定一個用ndarray表示的初始解; - complex_constraints:復雜約束，默認為None,表示沒有約束，具體定義同GA 的 complex_constraints; - complex_constraints_method:求解復雜約束的方法，默認為'loop',即如果解不滿足復雜約束，則再次隨機產生解，直到滿足約束,暫時不支持其他的求解方式。 RandomWalk 除了提供基本的random_walk函數之外，還提供了一個更加強大的improved_random_walk函數，后者的全局尋優能力要更強。

6. 求解帶復雜約束的目標函數

    上面所述的各種優化方法求解的都是變量僅有簡單邊界約束(形如$a \le x_i \le b$),下面介紹如何使用各種優化方法求解帶有復雜約束條件的目標函數。其實，求解方法也非常簡單，以GA為例，下面的例子即對Rosenbrock函數求解了帶有三個復雜約束條件的最優值: ```python from time import time from sopt.GA.GA import GA from sopt.util.functions import * from sopt.util.ga_config import * from sopt.util.constraints import * class TestGA: def __init__(self): self.func = Rosenbrock self.func_type = Rosenbrock_func_type self.variables_num = Rosenbrock_variables_num self.lower_bound = Rosenbrock_lower_bound self.upper_bound = Rosenbrock_upper_bound self.cross_rate = 0.8 self.mutation_rate = 0.1 self.generations = 300 self.population_size = 200 self.binary_code_length = 20 self.cross_rate_exp = 1 self.mutation_rate_exp = 1 self.code_type = code_type.real self.cross_code = False self.select_method = select_method.proportion self.rank_select_probs = None self.tournament_num = 2 self.cross_method = cross_method.uniform self.arithmetic_cross_alpha = 0.1 self.arithmetic_cross_exp = 1 self.mutation_method = mutation_method.uniform self.none_uniform_mutation_rate = 1 self.complex_constraints = [constraints1,constraints2,constraints3] self.complex_constraints_method = complex_constraints_method.penalty self.complex_constraints_C = 1e8 self.M = 1e8 self.GA = GA(**self.__dict__) def test(self): start_time = time() self.GA.run() print("GA costs %.4f seconds!" % (time()-start_time)) self.GA.save_plot() self.GA.show_result() if __name__ == '__main__': TestGA().test() ``` 運行結果如下: ``` GA costs 1.9957 seconds! -------------------- GA config is: -------------------- lower_bound:[-2.048, -2.048] cross_code:False complex_constraints_method:penalty mutation_method:uniform mutation_rate:0.1 mutation_rate_exp:1 cross_rate:0.8 upper_bound:[2.048, 2.048] arithmetic_cross_exp:1 variables_num:2 generations:300 tournament_num:2 select_method:proportion func_type:max complex_constraints_C:100000000.0 cross_method:uniform complex_constraints:[ , , ] func: none_uniform_mutation_rate:1 cross_rate_exp:1 code_type:real M:100000000.0 binary_code_length:20 global_generations_step:300 population_size:200 arithmetic_cross_alpha:0.1 -------------------- GA caculation result is: -------------------- global best target generation index/total generations:226/300 global best point:[ 1.7182846 -1.74504313] global best target:2207.2089435117955 ``` 圖6 GA求解帶有三個復雜約束的Rosenbrock函數 上面的constraints1,constraints2,constraints3是三個預定義的約束條件函數，其定義分別為:$constraints1:x_1^2 + x_2^2 - 6 \le 0$;$constraints2:x_1 + x_2 \le 0$;$constraints3:-2-x_1 - x_2 \le 0$,函數原型為: ```python def constraints1(x): x1 = x[0] x2 = x[1] return x1**2 + x2**2 -3 def constraints2(x): x1 = x[0] x2 = x[1] return x1+x2 def constraints3(x): x1 = x[0] x2 = x[1] return -2 -x1 -x2 ``` 其實觀察可以發現，上面的代碼和原始的GA實例代碼唯一的區別，就是其增加了`self.complex_constraints = [constraints1,constraints2,constraints3]`這樣一句，對於其他的優化方法，其都定義了`complex_constraints`和`complex_constraints_method`這兩個屬性，只要傳入相應的約束條件函數列表以及求解約束條件的方法就可以求解帶復雜約束的目標函數了。比如我們再用Random Walk求解和上面一樣的帶三個約束的Rosenbrock函數，代碼及運行結果如下: ```python from time import time from sopt.util.functions import * from sopt.util.constraints import * from sopt.util.random_walk_config import * from sopt.Optimizers.RandomWalk import RandomWalk class TestRandomWalk: def __init__(self): self.func = Rosenbrock self.func_type = Rosenbrock_func_type self.variables_num = Rosenbrock_variables_num self.lower_bound = Rosenbrock_lower_bound self.upper_bound = Rosenbrock_upper_bound self.generations = 100 self.init_step = 10 self.eps = 1e-4 self.vectors_num = 10 self.init_pos = None self.complex_constraints = [constraints1,constraints2,constraints3] self.complex_constraints_method = complex_constraints_method.loop self.RandomWalk = RandomWalk(**self.__dict__) def test(self): start_time = time() self.RandomWalk.random_walk() print("random walk costs %.4f seconds!" %(time() - start_time)) self.RandomWalk.save_plot() self.RandomWalk.show_result() if __name__ == '__main__': TestRandomWalk().test() ``` 運行結果: ``` Finish 1 random walk! Finish 2 random walk! Finish 3 random walk! Finish 4 random walk! Finish 5 random walk! Finish 6 random walk! Finish 7 random walk! Finish 8 random walk! Finish 9 random walk! Finish 10 random walk! Finish 11 random walk! Finish 12 random walk! Finish 13 random walk! Finish 14 random walk! Finish 15 random walk! Finish 16 random walk! Finish 17 random walk! random walk costs 0.1543 seconds! -------------------- random walk config is: -------------------- eps:0.0001 func_type:max lower_bound:[-2.048 -2.048] upper_bound:[2.048 2.048] init_step:10 vectors_num:10 func: variables_num:2 walk_nums:17 complex_constraints_method:loop generations:100 generations_nums:2191 complex_constraints:[ , , ] -------------------- random walk caculation result is: -------------------- global best generation index/total generations:2091/2191 global best point is: [-2.41416736 0.41430367] global best target is: 2942.6882849234585 ``` 圖7 Random Walk求解帶有三個復雜約束的Rosenbrock函數 可以發現Random Walk 求解得到的最優解要比GA好，而且運行時間更快，經過實驗發現，在所有的優化方法中，不論是求解帶復雜約束還是不帶復雜約束條件的目標函數，求解效果大體上排序是:Random Walk > PSO > GA > SA 。所以當你在求解具體問題時，不妨多試幾種優化方法，然后擇優選擇。

7. 基於梯度的系列優化方法

    上面所述的各種優化方法，比如GA,PSO,SA等都是基於隨機搜索的優化算法，其計算是不依賴於目標函數的具體形式，也不需要知道其梯度的，更加傳統的優化算法是基於梯度的算法，比如經典的梯度下降(上升)法(Gradient Descent)以及其一系列變種。下面就簡要介紹sopt中GD,Momentum,AdaGrad,RMSProp以及Adam的實現。關於這些基於梯度的優化算法的具體原理，可以參考我之前的一篇博文[深度學習中常用的優化方法](http://www.cnblogs.com/lyrichu/p/8940363.html)。另外需要注意，以下所述的基於梯度的各種優化算法，一般都是用在無約束優化問題里面的，如果是有約束的問題，請選擇上面其他的優化算法。下面是GradientDescent的一個使用實例: ```python from time import time from sopt.util.gradients_config import * from sopt.util.functions import * from sopt.Optimizers.Gradients import GradientDescent class TestGradientDescent: def __init__(self): self.func = quadratic50 self.func_type = quadratic50_func_type self.variables_num = quadratic50_variables_num self.init_variables = None self.lr = 1e-3 self.epochs = 5000 self.GradientDescent = GradientDescent(**self.__dict__) def test(self): start_time = time() self.GradientDescent.run() print("Gradient Descent costs %.4f seconds!" %(time()-start_time)) self.GradientDescent.save_plot() self.GradientDescent.show_result() if __name__ == '__main__': TestGradientDescent().test() ``` 運行結果為: ``` Gradient Descent costs 14.3231 seconds! -------------------- Gradient Descent config is: -------------------- func_type:min variables_num:50 func: epochs:5000 lr:0.001 -------------------- Gradient Descent caculation result is: -------------------- global best epoch/total epochs:4999/5000 global best point: [ 0.9999524 1.99991045 2.99984898 3.9998496 4.99977767 5.9997246 6.99967516 7.99964102 8.99958143 9.99951782 10.99947879 11.99944665 12.99942492 13.99935192 14.99932708 15.99925856 16.99923686 17.99921689 18.99911527 19.9991255 20.99908968 21.99899699 22.99899622 23.99887832 24.99883597 25.99885616 26.99881394 27.99869772 28.99869349 29.9986766 30.99861142 31.99851987 32.998556 33.99849351 34.99845985 35.99836731 36.99832444 37.99831792 38.99821067 39.99816567 40.99814951 41.99808199 42.99808161 43.99806655 44.99801207 45.99794449 46.99788003 47.99785468 48.99780825 49.99771656] global best target: 1.0000867498727912 ``` 圖7 GradientDescent 求解quadratic50 下面簡要說明以下GradientDescent,Momentum等類中的主要參數,像`func`,`variables_num`等含義已經解釋很多次了，不再贅述，這里主要介紹各類特有的一些參數。 1. GradientDescent類: - lr:學習率,默認是$10^{-3}$; - epochs:迭代次數，默認是100； 2. Momentum類: - lr:學習率,默認是$10^{-3}$; - beta:Momentum的$\beta$參數，默認是0.9; - epochs:迭代次數，默認是100； 3. AdaGrad類: - lr:學習率,默認是$10^{-3}$; - eps:極小的一個正數(防止分母為0),默認為$10^{-8}$; - epochs:迭代次數，默認是100； 4. RMSProp類: - lr:學習率,默認是$10^{-3}$; - beta:RMSProp的$\beta$參數，默認是0.9; - eps:極小的一個正數(防止分母為0),默認為$10^{-8}$; - epochs:迭代次數，默認是100； 5. Adam類: - lr:學習率,默認是$10^{-3}$; - beta1:Adam的$\beta_1$參數，默認是0.5; - beta2:Adam的$\beta_2$參數，默認是0.9; - eps:極小的一個正數(防止分母為0),默認為$10^{-8}$; - epochs:迭代次數，默認是100； 這里額外說一句,Adam原文作者給出的最優超參數為$\beta_1 = 0.9,\beta_2 = 0.999$,但是我在實際調試過程中發現,$\beta_1 = 0.5,\beta_2 = 0.9$才能達到比較好的效果，具體原因目前不是很清楚。如果想要學習各個優化函數類的具體用法，還可以使用`from sopt.test import *`導入預定義的示例測試，來運行觀察結果。比如下面的代碼就運行了PSO的一個示例測試: ```python from sopt.test import test_PSO test_PSO.TestPSO().test() ``` 結果: ``` PSO costs 3.4806 seconds! -------------------- PSO config is: -------------------- lower_bound:[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1] generations:200 vmin:-1 func: w_method:linear_decrease func_type:min population_size:100 w_start:0.9 complex_constraints_method:loop vmax:1 c2:1.49445 upper_bound:[50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50] variables_num:50 c1:1.49445 w:1 complex_constraints:None w_end:0.4 -------------------- PSO calculation result is: -------------------- global best generation index/total generations: 199/200 global best point: [ 1. 1. 2.94484328 4.13875216 5.00293498 6.13124759 6.99713025 7.92116383 8.87648843 10.02066994 11.0758768 12.02240279 13.01125368 13.98010373 14.98063168 15.97776149 17.11878537 18.00246112 18.14780887 20.00637617 21.00223704 22.00689373 23.14823218 24.0002456 24.98672157 25.99141686 27.02112321 28.01540506 29.05403155 30.07304888 31.00414822 32.00982867 32.99444884 33.9114213 34.96631157 36.22871824 37.0015616 37.98907918 39.01245751 40.1371835 41.0182043 42.07768102 42.87178292 43.93687997 45.05786395 46.03778693 47.07913415 50. 48.9964866 50. ] global best target: 6.95906097685 ```

附錄

    附錄會簡要說明上文提到的一些概念，待有空更新。。。

To do ...

1. GA的並行化計算 2. 小生境GA 3. GA,PSO等求解離散函數最優解 4. GA,PSO等實際應用舉例 5. 其他的一些最優化方法,比如牛頓法，擬牛頓法等 6. 其他