编写带有约束的用户定义问题（NLP）#
我们在这里展示如何编写一个具有单一目标、等式和不等式约束的非平凡用户定义问题（UDP）。 我们假设问题的数学公式如下：
这是Luksan, L., 和 Jan Vlcek在“无约束和等式约束优化的稀疏和部分可分离测试问题”（1999）中问题5.9的修改实例。修改在于最后两个约束，在本教程中，这些约束被视为不等式而非等式约束。
当前问题具有边界约束，4个等式约束，两个不等式约束（注意这些约束的不同形式）和一个目标函数，在优化问题的分类中，可以归类为非线性规划（NLP）问题。
暂时忽略适应度，UDP 让 pygmo 理解问题类型的基本结构将是：
>>> class my_constrained_udp:
... def get_bounds(self):
... return ([-5]*6,[5]*6)
... # Inequality Constraints
... def get_nic(self):
... return 2
... # Equality Constraints
... def get_nec(self):
... return 4
注意我们需要如何指定等式约束的数量和不等式约束的数量（因为pygmo默认假设两者都为0）。不需要指定目标的数量，因为pygmo默认假设为单目标优化。关于UDP规范的完整文档可以在pygmo.problem文档中找到。
我们仍然需要编写适应度函数，因为这是所有UDPs的强制方法（与get_bounds()一起）。使用不完整的UDP构建problem将会失败。在pygmo中，适应度包括目标和约束，按照描述的顺序[obj,ec,ic]。所有等式约束的形式为\(g(x) = 0\)，而不等式约束的形式为\(g(x) <= 0\)，如pygmo.problem.fitness()中所述。
>>> import math
>>> class my_constrained_udp:
... def fitness(self, x):
... obj = 0
... for i in range(3):
... obj += (x[2*i-2]-3)**2 / 1000. - (x[2*i-2]-x[2*i-1]) + math.exp(20.*(x[2*i - 2]-x[2*i-1]))
... ce1 = 4*(x[0]-x[1])**2+x[1]-x[2]**2+x[2]-x[3]**2
... ce2 = 8*x[1]*(x[1]**2-x[0])-2*(1-x[1])+4*(x[1]-x[2])**2+x[0]**2+x[2]-x[3]**2+x[3]-x[4]**2
... ce3 = 8*x[2]*(x[2]**2-x[1])-2*(1-x[2])+4*(x[2]-x[3])**2+x[1]**2-x[0]+x[3]-x[4]**2+x[0]**2+x[4]-x[5]**2
... ce4 = 8*x[3]*(x[3]**2-x[2])-2*(1-x[3])+4*(x[3]-x[4])**2+x[2]**2-x[1]+x[4]-x[5]**2+x[1]**2+x[5]-x[0]
... ci1 = 8*x[4]*(x[4]**2-x[3])-2*(1-x[4])+4*(x[4]-x[5])**2+x[3]**2-x[2]+x[5]+x[2]**2-x[1]
... ci2 = -(8*x[5] * (x[5]**2-x[4])-2*(1-x[5]) +x[4]**2-x[3]+x[3]**2 - x[4])
... return [obj, ce1,ce2,ce3,ce4,ci1,ci2]
... def get_bounds(self):
... return ([-5]*6,[5]*6)
... def get_nic(self):
... return 2
... def get_nec(self):
... return 4
... def gradient(self, x):
... return pg.estimate_gradient_h(lambda x: self.fitness(x), x)
为了检查上述UDP是否适合pygmo，我们尝试从中构造一个pygmo.problem并检查它：
>>> import pygmo as pg
>>> prob = pg.problem(my_constrained_udp())
>>> print(prob)
Problem name: ...
Global dimension: 6
Integer dimension: 0
Fitness dimension: 7
Number of objectives: 1
Equality constraints dimension: 4
Inequality constraints dimension: 2
Tolerances on constraints: [0, 0, 0, 0, 0, ... ]
Lower bounds: [-5, -5, -5, -5, -5, ... ]
Upper bounds: [5, 5, 5, 5, 5, ... ]
Has batch fitness evaluation: false
Has gradient: true
User implemented gradient sparsity: false
Expected gradients: 42
Has hessians: false
User implemented hessians sparsity: false
Fitness evaluations: 0
Gradient evaluations: 0
Thread safety: none
一切似乎都井然有序。尺寸符合我们的要求，检测到了梯度等。
解决您的受限UDP#
所以我们现在有一个带约束的UDP和一个数值梯度。让我们来解决它。可以部署许多不同的策略，我们这里只尝试两种：a) 使用增广拉格朗日方法 b) 使用单调盆地跳跃。考虑以下脚本：
>>> #METHOD A
>>> algo = pg.algorithm(uda = pg.nlopt('auglag'))
>>> algo.extract(pg.nlopt).local_optimizer = pg.nlopt('var2')
>>> algo.set_verbosity(200) # in this case this correspond to logs each 200 objevals
>>> pop = pg.population(prob = my_constrained_udp(), size = 1)
>>> pop.problem.c_tol = [1E-6] * 6
>>> pop = algo.evolve(pop)
objevals: objval: violated: viol. norm:
1 5.148e+30 6 203.761 i
201 1.27621 5 0.179321 i
401 1.71251 5 0.0550095 i
601 1.96643 5 0.0269182 i
801 2.21529 5 0.00340511 i
1001 2.25337 5 0.000478665 i
1201 2.25875 4 6.60584e-05 i
>>> print(pop.get_fevals())
3740
>>> print(pop.get_gevals())
3696
解决方案在这里找到，调用了3740次适应度函数（约1200次目标评估和约2500次约束评估）和3696次梯度（在我们的情况下，每次对应6次适应度评估，因为使用了pygmo.estimate_gradient_h()来数值估计梯度）。
>>> #METHOD B
>>> algo = pg.algorithm(uda = pg.mbh(pg.nlopt("slsqp"), stop = 20, perturb = .2))
>>> algo.set_verbosity(1) # in this case this correspond to logs each 1 call to slsqp
>>> pop = pg.population(prob = my_constrained_udp(), size = 1)
>>> pop.problem.c_tol = [1E-6] * 6
>>> pop = algo.evolve(pop)
Fevals: Best: Violated: Viol. Norm: Trial:
14 3.59547e+36 6 501.393 0 i
19 1.89716e+38 5 432.423 0 i
39 1.89716e+38 5 432.423 1 i
44 1.89716e+38 5 432.423 2 i
49 1.18971e+28 5 231.367 0 i
171 2.25966 0 0 0
176 2.25966 0 0 1
379 2.25966 0 0 2
384 2.25966 0 0 3
389 2.25966 0 0 4
682 2.25966 0 0 5
780 2.25966 0 0 6
1040 2.25966 0 0 7
1273 2.25966 0 0 8
1278 2.25966 0 0 9
1415 2.25966 0 0 10
1558 2.25966 0 0 11
1563 2.25966 0 0 12
1577 2.25966 0 0 13
1645 2.25966 0 0 14
1878 2.25966 0 0 15
2051 2.25966 0 0 16
2179 2.25966 0 0 17
2184 2.25966 0 0 18
2189 2.25966 0 0 19
2194 2.25966 0 0 20
>>> pop.problem.get_fevals()
2195
>>> pop.problem.get_gevals()
1320
这些运行中的两种策略都收敛到了值为2.25966的局部可行最小值。从不同的初始种群重复上述解决方案策略，有时也会得到1.60799的值。
如果不需要，请不要使用黑箱求解器#
我们通过讨论如何反对科学界部分常见的（不良）实践来结束本教程，使用适当的局部优化算法总是更可取的，而启发式方法应仅在需要的情况下使用，否则它们只是一个糟糕的想法。让我们考虑一下在CEC2006竞赛期间使用的约束优化问题套件。在pygmo中，我们在类pygmo.cec2006中实现了这样的UDP。由于竞赛旨在进行启发式优化，因此该类没有实现梯度。因此，我们快速编写了一个元UDP，为没有梯度的UDP添加数值梯度：
>>> class add_gradient:
... def __init__(self, prob):
... self.prob = pg.problem(prob)
... def fitness(self, x):
... return self.prob.fitness(x)
... def get_bounds(self):
... return self.prob.get_bounds()
... def get_nec(self):
... return self.prob.get_nec()
... def get_nic(self):
... return self.prob.get_nic()
... def get_nobj(self):
... return self.prob.get_nobj()
... def gradient(self, x):
... return pg.estimate_gradient(lambda x: self.fitness(x), x) # we here use the low precision gradient
超级酷，我知道。这样的元UDP非常有用，因为它现在允许在CEC2006问题实例上调用例如SQP方法。
这里只考虑一个案例：问题g07。你可以通过研究其余案例中发生的情况来完成本教程。
>>> # Start adding the numerical gradient (low-precision) to the UDP
>>> prob = pg.problem(add_gradient(pg.cec2006(prob_id = 7)))
>>> # Define a solution strategy (SQP method)
>>> algo = pg.algorithm(uda = pg.mbh(pg.nlopt("slsqp"), 20, .2))
>>> pop = pg.population(prob, 1)
>>> # Solve the problem
>>> pop = algo.evolve(pop)
>>> # Collect information
>>> print(pop.champion_f[0])
24.306209068189599
>>> pop.problem.get_fevals()
1155
>>> pop.problem.get_gevals()
1022
因此，解决问题的总评估次数（精度为1e-8）为1155 + 1022 * 2 = 3599。 为了比较，作为一个例子，我们检查了启发式方法可能提供的表格，该表格出现在：
黄, 薇琪·凌, A. 凯·秦, 和 波努图拉伊·N·苏甘坦. “用于约束实数参数优化的自适应差分进化算法.” 进化计算, 2006. CEC 2006. IEEE 大会. IEEE, 2006.
发现经过5000次适应度评估后，本文介绍的启发式方法仍未解决这个特定问题。