编写一个简单的用户定义问题#

虽然pagmo提供了许多UDP（参见问题列表）来帮助您测试自己的优化策略或用户定义的算法，但编写自己的UDP是pygmo使用的核心。在本教程中，我们将展示如何编写UDP。请记住，UDP是可以用来构建问题的类，而算法可以解决这些问题。请参阅使用算法类和使用问题类以了解这些核心类的使用。

我们鼓励用户阅读类问题的文档，以获取可以在UDP中实现或必须实现的详细方法列表。为了简单起见，我们考虑最小化二维球体函数的简单问题。

\[\begin{split}\begin{array}{ll} \mbox{minimize: } & x_1^2+x_2^2 \\ \mbox{subject to:} & -1 \le x_i \le 1, i = 1..2 \end{array}\end{split}\]

注意

在pagmo中，总是假设为最小化，如果你需要最大化某个目标函数，只需在该目标前加上一个负号。

>>> class sphere_function:
...     def fitness(self, x):
...         return [sum(x*x)]
...
...     def get_bounds(self):
...         return ([-1,-1],[1,1])

您必须在类中实现的两个必需方法是 fitness(self, x) 和 get_bounds(self)。问题维度将通过第二个方法的返回值推断，而决策向量的实际适应度将通过调用第一个方法并将 x 作为 NumPy 数组传递来计算。重要的是要记住 x 是一个 NumPy 数组，因此 NumPy 数组算术适用于 fitness() 的主体。还要注意，定义 UDP 时我们不需要从其他 pygmo 相关类继承。由于在这种情况下我们没有定义任何其他方法，pygmo 将假设为单目标、无约束、无梯度等…

现在让我们从新的UDP构建一个问题。

>>> import pygmo as pg
>>> prob = pg.problem(sphere_function())

这很简单！要检查pygmo从我们的UDP中检测到的问题类型，我们可以在屏幕上打印：

>>> print(prob) 
Problem name: <class 'sphere_function'>
    C++ class name: ...

    Global dimension:                       2
    Integer dimension:                      0
    Fitness dimension:                      1
    Number of objectives:                   1
    Equality constraints dimension:         0
    Inequality constraints dimension:       0
    Lower bounds: [-1, -1]
    Upper bounds: [1, 1]
    Has batch fitness evaluation: false

    Has gradient: false
    User implemented gradient sparsity: false
    Has hessians: false
    User implemented hessians sparsity: false

    Fitness evaluations: 0

    Thread safety: none

现在让我们增加一些（轻微的）复杂性。我们希望我们的UDP具有可扩展性：

\[\begin{split}\begin{array}{ll} \mbox{minimize: } & \sum_i x_i^2 \\ \mbox{subject to:} & -1 \le x_i \le 1, i = 1..n \end{array}\end{split}\]

并且有一个人类可读的名称。

>>> class sphere_function:
...     def __init__(self, dim):
...         self.dim = dim
...
...     def fitness(self, x):
...         return [sum(x*x)]
...
...     def get_bounds(self):
...         return ([-1] * self.dim, [1] * self.dim)
...
...     def get_name(self):
...         return "Sphere Function"
...
...     def get_extra_info(self):
...         return "\tDimensions: " + str(self.dim)
>>> prob = pg.problem(sphere_function(3))
>>> print(prob) 
Problem name: Sphere Function
    C++ class name: ...

    Global dimension:                       3
    Integer dimension:                      0
    Fitness dimension:                      1
    Number of objectives:                   1
    Equality constraints dimension:         0
    Inequality constraints dimension:       0
    Lower bounds: [-1, -1, -1]
    Upper bounds: [1, 1, 1]
    Has batch fitness evaluation: false

    Has gradient: false
    User implemented gradient sparsity: false
    Has hessians: false
    User implemented hessians sparsity: false

    Fitness evaluations: 0

    Thread safety: none

Extra info:
    Dimensions: 3

好吧，那很简单，但现在有一个问题需要解决…

>>> algo = pg.algorithm(pg.bee_colony(gen = 20, limit = 20))
>>> pop = pg.population(prob,10)
>>> pop = algo.evolve(pop)
>>> print(pop.champion_f) 
[  3.75822114e-06]

哇，那些蜜蜂！！

可能的陷阱#

嗯，这很好，因为它像魔法一样有效。但UDP也可能是一个相当复杂的类，它有可能以某种方式格式错误。让我们看看一些常见的错误。

>>> class sphere_function:
...     def fitness(self, x):
...         return [sum(x*x)]
...
>>> pg.problem(sphere_function()) 
NotImplementedError                       Traceback (most recent call last)
...
NotImplementedError: the mandatory 'get_bounds()' method has not been detected in the user-defined Python problem
'<sphere_function object at 0x1108cad68>' of type '<class 'sphere_function'>': the method is either not present or not callable

哎呀，我忘记实现两个必需方法中的一个了。在这种情况下，无法构造一个问题，当我们尝试时，会得到一个非常有用的错误信息。

在其他情况下，尽管UDP仍然是格式错误的，问题的构建将会成功，问题只有在调用格式错误的方法时才会显现：

>>> class sphere_function:
...     def fitness(self, x):
...         return sum(x*x)
...
...     def get_bounds(self):
...         return ([-1,-1],[1,1])
>>> prob = pg.problem(sphere_function())
>>> prob.fitness([1,2]) 
AttributeError                            Traceback (most recent call last)
...
AttributeError: 'numpy.float64' object has no attribute '__iter__'

在这种情况下，问题是fitness()方法返回的是一个标量而不是一个类似数组的对象（记住pygmo也能够解决多目标和约束问题，因此适应度值通常必须是一个向量）。pygmo在第一次调用fitness()方法时会抱怨错误的返回类型。

计算速度的注意事项#

编写UDP的最高效方式是使用C++编写并将其暴露给Python。pygmo中包含的大多数UDP（参见问题列表）都是这样实现的。然而，在编写自己的UDP时，通常更快且更少痛苦的方式是直接在Python中编写，如本教程所示。这对理想情况有什么影响呢？让我们看看，在一台测试机器上，一个简单的例子：可扩展的Rosenbrock函数：

\[\begin{split}\begin{array}{ll} \mbox{minimize: } & \sum_{i=1}^{N-1} 100 (x_{i+1} - x_i^2 )^2 + (1-x_i)^2 \\ \mbox{subject to:} & -5 \le x_i \le 10, i = 1..N \end{array}\end{split}\]

在pygmo中可以快速写成：

>>> import numpy as np
>>> class py_rosenbrock:
...     def __init__(self,dim):
...         self.dim = dim
...     def fitness(self,x):
...         retval = np.zeros((1,))
...         for i in range(len(x) - 1):
...             retval[0] += 100.*(x[i + 1]-x[i]**2)**2+(1.-x[i])**2
...         return retval
...     def get_bounds(self):
...         return (np.full((self.dim,),-5.),np.full((self.dim,),10.))

我们现在快速且粗略地分析一个高维Rosenbrock实例的实例化：2000个变量！！

>>> prob_python = pg.problem(py_rosenbrock(2000))
>>> prob_cpp = pg.problem(pg.rosenbrock(2000))
>>> dummy_x = np.full((2000,), 1.)
>>> import time
>>> start_time = time.time(); [prob_python.fitness(dummy_x) for i in range(1000)]; print(time.time() - start_time) 
2.3352...
>>> start_time = time.time(); [prob_cpp.fitness(dummy_x) for i in range(1000)]; print(time.time() - start_time) 
0.0114226...

等一下……真的吗？Python比C++慢两个数量级？别慌。这是一个非常大的问题，而且那个for循环在Python中不会超级优化。让我们看看在这些情况下我们是否能做得更好……让我们使用numba中的jit装饰器将我们的适应度方法编译成C代码。

>>> from numba import jit
>>> class jit_rosenbrock:
...     def __init__(self,dim):
...         self.dim = dim
...     @jit
...     def fitness(self,x):
...         retval = np.zeros((1,))
...         for i in range(len(x) - 1):
...             retval[0] += 100.*(x[i + 1]-x[i]**2)**2+(1.-x[i])**2
...         return retval
...     def get_bounds(self):
...         return (np.full((self.dim,),-5.),np.full((self.dim,),10.))
>>> prob_jit = pg.problem(jit_rosenbrock(2000))
>>> start_time = time.time(); [prob_jit.fitness(dummy_x) for i in range(1000)]; print(time.time() - start_time) 
0.03771...

通过更多的努力，我们可以进一步提高性能：

>>> from numba import jit, float64
>>> class jit_rosenbrock2:
...      def fitness(self,x):
...          return jit_rosenbrock2._fitness(x)
...      @jit(float64[:](float64[:]),nopython=True)
...      def _fitness(x):
...          retval = np.zeros((1,))
...          for i in range(len(x) - 1):
...              tmp1 = (x[i + 1]-x[i]*x[i])
...              tmp2 = (1.-x[i])
...              retval[0] += 100.*tmp1*tmp1+tmp2*tmp2
...          return retval
...      def get_bounds(self):
...          return (np.full((self.dim,),-5.),np.full((self.dim,),10.))
...      def __init__(self,dim):
...          self.dim = dim
>>> prob_jit2 = pg.problem(jit_rosenbrock2(2000))
>>> start_time = time.time(); [prob_jit2.fitness(dummy_x) for i in range(1000)]; print(time.time() - start_time) 
0.01687...

好多了，对吧？

注意

有关使用Numba加速Python代码的更多信息，请参阅Numba文档页面。特别是，请注意，此用途仅支持NumPy和Python语言的一部分。

编写一个简单的用户定义问题

目录

编写一个简单的用户定义问题#

可能的陷阱#

计算速度的注意事项#