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$,其中$02. 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. 其他