Mathematica中的约束最优化问题/遗传算法/退火算法

引言
约束非线性最优化的数值算法大致可以分为 基于梯度的方法 和 直接搜索方法. 基于梯度的方法使用一阶导数 (梯度) 或二阶导数（黑塞（Hessians）矩阵）. 例如，序列二次规划方法(SQP)、增广拉格朗日方法和（非线性）内点法. 直接搜索方法不使用导数信息. 例如，Nelder-Mead、遗传算法和差分进化，以及模拟退火. 直接搜索方法往往收敛速度较慢，但可以对函数和约束条件中噪声存在的限制更为宽松.

通常情况下，算法只建立问题的局部模型. 此外，许多这种算法总是寻求目标函数，或者由目标函数和约束条件合成的一个优值函数的某种减小，来确保迭代过程的收敛. 如果收敛的话，这种算法只找到局部最优解，因此被称为 局部最优化算法. 在 Wolfram 语言里，局部最优化问题可以用 FindMinimum 求解.

然而，全局优化算法 通常通过既允许减小也允许增大目标函数/优值函数来试图求全局最优解. 通常这种算法的计算代价更高. 全局最优化问题可以用 Minimize 精确求解或者用 NMinimize 数值化地求解.

这里我们求解一个非线性规划问题，

利用 Minimize，能够产生一个精确的解

In[1]:= Minimize[{x - y, -3 x^2 + 2 x y - y^2 >= -1}, {x, y}] https://wolfram.com/xid/0gmpon34wjytlihky4i-eu6i
Out[1]=
下面我们数值化地解同样的问题. NMinimize 给我们返回一个具有机器精度的解：
In[2]:= NMinimize[{x - y, -3 x^2 + 2 x y - y^2 >= -1}, {x, y}] https://wolfram.com/xid/0gmpon34wjytlihky4i-w8xkl
Out[2]=
FindMinimum 数值化地求一个局部最小值. 在这个例子中找到的局部最小值也是一个全局最小值.
In[3]:= FindMinimum[{x - y, -3 x^2 + 2 x y - y^2 >= -1}, {x, y}] https://wolfram.com/xid/0gmpon34wjytlihky4i-eveio
Out[3]=
NMinimize 函数
NMinimize 和 NMaximize 执行好几种算法用以求约束全局最优解. 这些方法具有足够的灵活性，以处理不可微或连续，并且不容易被局部最优解捕捉的函数.

寻找一个全局最优解的问题可具有任意难度，即使在没有约束的情况下，所以我们使用的方法可能会失败. 在通常情况下，多次使用不同的初始条件最优化，并采用其中最好的结果来，对于最优化该函数是有益的.

下面求 的最大值：
In[46]:= NMaximize[Sin[x + y] - x^2 - y^2, {x, y}] https://wolfram.com/xid/0gmpon34wjytlihky4i-cwooy7
Out[46]=
下面求 的最小值，约束条件为 和 ：
In[47]:= NMinimize[{x^2 + (y - .5)^2, y >= 0 && y >= x + 1}, {x, y}] https://wolfram.com/xid/0gmpon34wjytlihky4i-h7nrp3
Out[47]=
NMinimize 和 NMaximize 的约束条件可以是一个列表，或者是等式、不等式和定义域的划定的逻辑组合. 等式和不等式可能是非线性的. 由于近似数字处理上的局限，任何强不等式将被转换为弱不等式. 划定一个变量的定义域时用 Element ，例如，Element[x,Integers] 或者 x∈Integers. 变量必须或者是整数或者是实数，并且除非特别说明，否则将被假定为实数. 当点离开可行域时，约束条件通过增加罚函数进行强化.

约束条件可以包含如 And、Or 等等的逻辑运算符.
In[3]:= NMinimize[{x^2 + y^2, x >= 1 || y >= 2}, {x, y}] https://wolfram.com/xid/0gmpon34wjytlihky4i-b01n24
Out[3]=
这里 被限制为一个整数：
In[4]:= NMinimize[{(x - 1/3)^2 + (y - 1/3)^2, x [Element] Integers}, {x, y}] https://wolfram.com/xid/0gmpon34wjytlihky4i-hg9mh4
Out[4]=
为了让 NMinimize 工作，需要从一个长方形的初始区域内开始. 这和其他数值方法中需要提供一个初始点或几个初始点类似. 初始区域是通过给每个变量一个有限的上界和下界来确定的. 这通过在约束条件中包括 ，或者在变量中包括 {x,a,b} 完成. 如果两者都被提供，则变量中的边界用于初始区域，而约束条件就作为约束条件来使用. 如果对于一个变量 x 没有指定初始区域，我们将采用默认初始区域 . 不同变量可以有用不同方式定义的初始区域.

下面的初始区域由变量提供. 该问题是无约束的.
In[5]:= NMinimize[x^2, {{x, 3, 4}}] https://wolfram.com/xid/0gmpon34wjytlihky4i-fpmco8
Out[5]=
下面的初始区域由约束条件提供：
In[6]:= NMinimize[{x^2, 3 <= x <= 4}, {x}] https://wolfram.com/xid/0gmpon34wjytlihky4i-cxvs1g
Out[6]=
这里， 的初始区域由约束条件提供， 的初始区域由变量提供，而 的初始区域采用默认值. 该问题在 和 上不受约束，但是在 上则不然：
In[7]:= NMinimize[{x^2 + y^2 + z^2, 3 <= x <= 4}, {x, {y, 2, 5}, z}] https://wolfram.com/xid/0gmpon34wjytlihky4i-c0gzbf
Out[7]=
多项式 在 有全局极小值. NelderMead 找到极小值之一：
In[48]:= NMinimize[4 x^4 - 4 x^2 + 1, x, Method -> “NelderMead”] https://wolfram.com/xid/0gmpon34wjytlihky4i-n7i88i
Out[48]=
其他的最小值可以使用一个不同的 RandomSeed 来找到：
In[50]:= NMinimize[4 x^4 - 4 x^2 + 1, x,
Method -> {“NelderMead”, “RandomSeed” -> 111}] https://wolfram.com/xid/0gmpon34wjytlihky4i-o0tu53
Out[50]=
NMinimize 和 NMaximize 有好几个优化方法可采用，包括: Automatic、“DifferentialEvolution”、“NelderMead”、“RandomSearch” 和 “SimulatedAnnealing”. 优化方法由 Method 选项控制，选项可以取作为一个字符串的方法，也可以取一个列表，该列表的第一个元素是作为一个字符串的方法而其他元素都是基于具体方法的选项. 所有基于具体方法的选项，左边也必须作为字符串给出.

下面的函数有大量的局部极小值.
In[34]:= Clear[f];
f = E^Sin[50 x] + Sin[60 E^y] + Sin[70 Sin[x]] + Sin[Sin[80 y]] -
Sin[10 (x + y)] + 1/4 (x^2 + y^2);
Plot3D[f, {x, -1, 1}, {y, -1, 1}, PlotPoints -> 50, Mesh -> False] https://wolfram.com/xid/0gmpon34wjytlihky4i-bbktwe
Out[36]=
采用 RandomSearch 来寻找极小值：
In[54]:= NMinimize[f, {x, y}, Method -> “RandomSearch”] https://wolfram.com/xid/0gmpon34wjytlihky4i-fe5jvw
Out[54]=
采用具有更多初始值的 RandomSearch 来寻找全局极小值：
In[55]:= NMinimize[f, {x, y},
Method -> {“RandomSearch”, “SearchPoints” -> 250}] https://wolfram.com/xid/0gmpon34wjytlihky4i-doh666
Out[55]=
使用默认的方法，NMinimize 根据问题的类型选择使用的方法. 如果目标函数和约束条件是线性的，则采用LinearProgramming. 如果有整数变量，或者如果目标函数的头部不是一个数值函数的话，则采用差分进化算法. 对于所有其他类型的问题，则采用Nelder-Mead，但是如果Nelder-Mead 运行不佳的话，就切换到差分进化算法.

因为 NMinimize 所采用的方法不可能改善每次迭代的结果，所以我们只在数次迭代后检查收敛度.

约束全局最优化的数值算法
Nelder–Mead
Nelder-Mead方法是一个直接搜索方法. 对于一个具有 个变量的函数，该算法保持一个 个点的点集，这些点形成 维空间的多面体的顶点. 这种方法通常被称为 “单纯” 方法，它不应与众所周知的线性规划的单纯型法混为一谈.

在每次迭代中， 个点 形成一个多面体. 点的排列顺序使得 然后产生一个新的点以取代最差点

设 为包含最好的 个点的多面体的中心，. 一个试验点 由穿过中心的最差点反射产生， ，其中 是一个参数.

如果新的点 既不是一个新的最差点也不是一个新的最佳点，, 则由 取代 .

如果新的点 优于最佳点，，则反射非常成功并可以进一步执行到 ，其中 是扩大多面体的一个参数. 如果扩大成功，，则由 取代 ；否则扩大失败，并由 取代 .

如果新的点 劣于第二差点，，则认为多面体过于庞大，需要收缩. 由

定义一个新的尝试点，其中 是一个参数. 如果 ，收缩是成功的，且由 取代 . 否则就进行进一步的收缩.

如果在新的和旧的多面体里的最佳函数值的差，以及在新的最佳点和旧的最佳点的距离小于AccuracyGoal 和 PrecisionGoal 所提供的容差， 这个过程则被认为是收敛的.

严格地说，Nelder-Mead不是一个真正的全局最优化算法，然而在实践中对于没有许多局部极小值的问题，往往有不错的工作效率.

选项名
默认值
“ContractRatio” 0.5 用于收缩的比例
“ExpandRatio” 2.0 用于扩大的比例
“InitialPoints” Automatic 初始点集合
“PenaltyFunction” Automatic 应用于约束条件以惩罚无效点的函数
“PostProcess” Automatic 是否采用局部搜索方法进行后处理
“RandomSeed” 0 随机数生成器的初始值
“ReflectRatio” 1.0 用于反射的比例
“ShrinkRatio” 0.5 用于收缩的比例
“Tolerance” 0.001 接受对约束违反程度的宽限
NelderMead 具体选项.

用 NelderMead 对单位圆盘里的函数最小化.
In[82]:= https://wolfram.com/xid/0gmpon34wjytlihky4i-bznj2y
Out[82]=
这里是一个具有若干全部在不同深度的局部极小值的函数：
In[37]:= Clear[a, f]; a =
Reverse /@
Distribute[{{-32, -16, 0, 16, 32}, {-32, -16, 0, 16, 32}},
List]; f =
1/(0.002 +
Plus @@ MapIndexed[1/(#2[[1]] + Plus @@ (({x, y} - #1)^6)) &,
a]); Plot3D[f, {x, -50, 50}, {y, -50, 50}, Mesh -> None,
PlotPoints -> 25] https://wolfram.com/xid/0gmpon34wjytlihky4i-brkltd
Out[37]=
用默认参数，NelderMead 很容易陷入一个局部极小值里：
In[116]:= Do[Print[NMinimize[f, {{x, -50, 50}, {y, -50, 50}},
Method -> {“NelderMead”, “RandomSeed” -> i}]], {i, 5}] https://wolfram.com/xid/0gmpon34wjytlihky4i-02wy7

通过使用更积极并且更不可能使单纯型变小的设置，其结果便更好一些.
In[117]:= Do[Print[NMinimize[f, {{x, -50, 50}, {y, -50, 50}},
Method -> {“NelderMead”, “ShrinkRatio” -> 0.95,
“ContractRatio” -> 0.95, “ReflectRatio” -> 2,
“RandomSeed” -> i}]], {i, 5}] https://wolfram.com/xid/0gmpon34wjytlihky4i-bchagl

差分进化
差分进化是一个简单的随机函数最小化算法.

该算法保持了一个 个点的集合，，其中通常 ，这里 是变量的数目.

在算法的每次迭代过程中，都产生一个新的 个点的集合. 第 个新的点通过从旧的集合选择三个随机点 、 和 ，并且形成 来产生，其中 是一个实数尺度因子. 然后从 和 构建一个新的点 ，其中以概率 从 选择第 个坐标， 其它则从 选择坐标. 如果 ，那么就由 取代集合中的 . 概率 通过"CrossProbability" 选项控制.

如果在新的和旧的集合里最佳函数值的差，以及在新的最佳点和旧的最佳点的距离小于AccuracyGoal 和 PrecisionGoal 所提供的容差， 这个过程则被认为是收敛的.

差分进化方法在计算上是昂贵的，但是相对来说也较为稳健，而且对于有较多局部极小值的问题来说工作性能较好.

选项名
默认值
“CrossProbability” 0.5 从 xi 获取一个基因的概率
“InitialPoints” Automatic 初始点的集合
“PenaltyFunction” Automatic 用于约束条件以惩罚无效点的函数
“PostProcess” Automatic 是否用局部搜索方法进行后处理
“RandomSeed” 0 随机数生成器的初始值
“ScalingFactor” 0.6 应用于创建一个相伴向量的差异向量的尺度
“SearchPoints” Automatic 用于进化的群体大小
“Tolerance” 0.001 接受对约束违反程度的宽限
DifferentialEvolution 具体选项.

这里通过利用 DifferentialEvolution 将单位圆盘里的函数最小化：
In[125]:= NMinimize[{100 (y - x²⁾2 + (1 - x)^2, x^2 + y^2 <= 1}, {x, y},
Method -> “DifferentialEvolution”] https://wolfram.com/xid/0gmpon34wjytlihky4i-lhcmm3
Out[125]=
下面的约束优化问题有一个全局最小值 7.667180068813135`：
In[126]:= Clear[f, c, v, x1, x2, y1, y2, y3] https://wolfram.com/xid/0gmpon34wjytlihky4i-bz4soq
In[127]:= f = 2 x1 + 3 x2 + 3 y1/2 + 2 y2 - y3/2;
c = {x1^2 + y1 == 5/4, x2^(3/2) + 3 y2/2 == 3, x1 + y1 <= 8/5,
4 x2/3 + y2 <= 3, y3 <= y1 + y2, 0 <= x1 <= 10, 0 <= x2 <= 10,
0 <= y1 <= 1, 0 <= y2 <= 1,
0 <= y3 <= 1, {y1, y2, y3} [Element] Integers
};
v = {x1, x2, y1, y2, y3}; https://wolfram.com/xid/0gmpon34wjytlihky4i-ecb2wp
采用 DifferentialEvolution 的默认设置，产生了一个不令人满意的解：
In[130]:= NMinimize[{f, c}, v, Method -> “DifferentialEvolution”] https://wolfram.com/xid/0gmpon34wjytlihky4i-bjtn42
Out[130]=
通过调整 ScalingFactor，结果要好得多. 在这种情况下，针对整数变量，增大了的 ScalingFactor 给 DifferentialEvolution 提供了更好的灵活性.
In[131]:= NMinimize[{f, c}, v,
Method -> {“DifferentialEvolution”, “ScalingFactor” -> 1}] https://wolfram.com/xid/0gmpon34wjytlihky4i-bzp8uj
Out[131]=
模拟退火算法
模拟退火算法是一个简单的随机函数最小化算法. 这是从退火的物理过程得到启发的. 退火时，金属物体被加热到很高的温度并且慢慢冷却. 这个过程使金属的原子结构最终形成于一个较低的能量状态，从而成为更强硬的金属. 利用优化术语，退火使结构脱离局部最小值，探索并达到一个更好的，可望是全局的的最小值.

在每次迭代里，一个新的点，，在当前点 的附近产生. 临近区域的半径随每次迭代下降. 到目前为止，找到的最佳点，，也被追踪.

如果 ，则 取代 和 . 否则， 以概率 取代 . 这里 是由 BoltzmannExponent 定义的函数， 是当前的迭代， 是目标函数值的改变量，而 是上一次迭代中的目标函数值. 的默认函数是 .

正如 RandomSearch 方法，SimulatedAnnealing 使用多个初始点，并且从每个初始点出发找到一个最优值.

由选项 SearchPoints 给出的初始点数目的默认值为 ，其中 是变量数目.

对于每个初始点，如此反复执行直至达到最大迭代次数，该方法要么收敛到一个点，要么在由 LevelIterations 所提供的迭代次数里连续停留在同一点.

选项名
默认值
“BoltzmannExponent” Automatic 概率函数的指数
“InitialPoints” Automatic 初始点集
“LevelIterations” 50 停留在某一给定点的最大迭代次数
“PenaltyFunction” Automatic 用于约束条件以惩罚无效点的函数
“PerturbationScale” 1.0 随机跳转的尺度
“PostProcess” Automatic 是否利用局部搜索方法进行后处理
“RandomSeed” 0 随机数字生成器的初始值
“SearchPoints” Automatic 初始点的数目
“Tolerance” 0.001 接受对约束违反程度的宽限
SimulatedAnnealing 具体选项.

用 SimulatedAnnealing 对一个具有两个变量的函数最小化：
In[62]:= NMinimize[{100 (y - x²⁾2 + (1 - x)^2, -2.084 <= x <=
2.084 && -2.084 <= y <= 2.084}, {x, y},
Method -> “SimulatedAnnealing”] https://wolfram.com/xid/0gmpon34wjytlihky4i-klyhy
Out[62]=
下面是一个有多个局部最小值的函数：
In[38]:= Clear[f]
f[x_, y_] :=
20 Sin[[Pi]/2 (x - 2 [Pi])] +
20 Sin[[Pi]/2 (y - 2 [Pi])] + (x - 2 [Pi])^2 + (y - 2 [Pi])^2;
Plot3D[f[x, y], {x, 0, 10}, {y, 0, 10}] https://wolfram.com/xid/0gmpon34wjytlihky4i-ls9cpy
Out[40]=
默认情况下，SimulatedAnnealing 的步长没有大到足够摆脱局部极小值.
In[68]:= NMinimize[f[x, y], {x, y}, Method -> “SimulatedAnnealing”] https://wolfram.com/xid/0gmpon34wjytlihky4i-ffh8v2
Out[68]=
通过增加 PerturbationScale，采取的步长更大，产生的解更好.
In[69]:= https://wolfram.com/xid/0gmpon34wjytlihky4i-e9cyz9
Out[69]=
这里 BoltzmannExponent 被设置为使用一个能产生快速收敛的指数冷却函数.（请注意，修正的 PerturbationScale 仍然在使用.）
In[70]:= NMinimize[f[x, y], {x, y},
Method -> {“SimulatedAnnealing”, “PerturbationScale” -> 3,
“BoltzmannExponent” -> Function[{i, df, f0}, -df/(Exp[i/10])]}] https://wolfram.com/xid/0gmpon34wjytlihky4i-ef6al7
Out[70]=
随机搜索
随机搜索算法的工作方式是先产生一族随机初始点，然后从每个初始点用一个局部优化方法收敛到一个局部极小值. 最佳局部极小值被选为问题的解.

可能的局部搜索方法有 Automatic 和 “InteriorPoint”. 默认的方法是 Automatic，它使用 FindMinimum 把非约束的方法用于一个加了惩罚项的系统，而惩罚项是因约束而添加的. 当 Method 被设置为 “InteriorPoint” 时，则采用非线性内点方法.

由选项 SearchPoints 给定的初始点数目的默认值是 ，其中 是变量的数目.

RandomSearch 的收敛取决于每个初始点的局部方法的收敛度.

RandomSearch 运行很快，但其规模和搜索空间维度的匹配不是很好. 它也受到许多和 FindMinimum 相同的局限. 它不适合离散问题和那些导数或者secants不能提供多少有用信息的问题.

选项名
默认值
“InitialPoints” Automatic 初始点的集合
“Method” Automatic 用于最小化的方法
“PenaltyFunction” Automatic 应用于约束条件以惩罚无效点的函数
“PostProcess” Automatic 是否采用局部搜索方法进行后处理
“RandomSeed” 0 随机数字生成器的初始值
“SearchPoints” Automatic 用于启动局部搜索的点的数目
“Tolerance” 0.001 接受对约束违反程度的宽限
RandomSearch 具体选项.

这里利用 RandomSearch 将单位圆盘里的函数最小化：
In[71]:= NMinimize[{100 (y - x²⁾2 + (1 - x)^2, x^2 + y^2 <= 1}, {x, y},
Method -> “RandomSearch”] https://wolfram.com/xid/0gmpon34wjytlihky4i-diky5u
Out[71]=
这里是一个有好几个深度各不相同的局部极小值，通常难以优化的函数：
In[41]:= Clear[a, f]; a =
Reverse /@
Distribute[{{-32, -16, 0, 16, 32}, {-32, -16, 0, 16, 32}},
List]; f =
1/(0.002 +
Plus @@ MapIndexed[1/(#2[[1]] + Plus @@ (({x, y} - #1)^6)) &,
a]); Plot3D[f, {x, -50, 50}, {y, -50, 50}, Mesh -> None,
NormalsFunction -> “Weighted”, PlotPoints -> 50] https://wolfram.com/xid/0gmpon34wjytlihky4i-h70jdd
Out[41]=
采用 SearchPoints 的默认数目，有时并不能找到最小值：
In[73]:= Do[Print[NMinimize[f, {{x, -50, 50}, {y, -50, 50}},
Method -> {“RandomSearch”, “RandomSeed” -> i}]], {i, 5}] https://wolfram.com/xid/0gmpon34wjytlihky4i-d1rayw

用更多 SearchPoints 则产生较好好的结果：
In[74]:= Do[Print[NMinimize[f, {{x, -50, 50}, {y, -50, 50}},
Method -> {“RandomSearch”, “SearchPoints” -> 100,
“RandomSeed” -> i}]], {i, 5}] https://wolfram.com/xid/0gmpon34wjytlihky4i-c06fs

这里在一个网格上产生点作为初始点：
In[75]:= NMinimize[f, {{x, -50, 50}, {y, -50, 50}},
Method -> {“RandomSearch”,
“InitialPoints” ->
Flatten[Table[{i, j}, {i, -45, 45, 5}, {j, -45, 45, 5}], 1]}] https://wolfram.com/xid/0gmpon34wjytlihky4i-ba1r5o
Out[75]=
这里采用非线性内点法来寻找平方和的最小值：
In[76]:= n = 10;
f = Sum[(x[i] - Sin[i])^2, {i, 1, n}];
c = Table[-0.5 < x[i] < 0.5, {i, n}];
v = Array[x, n];
Timing[NMinimize[{f, c}, v,
Method -> {“RandomSearch”, Method -> “InteriorPoint”}]] https://wolfram.com/xid/0gmpon34wjytlihky4i-k6uuhk
Out[80]=
对于某些类型的问题，限制 SearchPoints 的数目可以更快地得到结果，并且不影响解的质量：
In[81]:= Timing[NMinimize[{f, c}, v,
Method -> {“RandomSearch”, Method -> “InteriorPoint”,
“SearchPoints” -> 1}]] https://wolfram.com/xid/0gmpon34wjytlihky4i-gh7ueg
Out[81]=