公司新闻

(六)最优化建模与算法之约束优化

约束优化问题在实际应用极其广泛，因为在实际问题中往往都存在诸多现实约束。约束优化考虑如下问题： $\\mathop{\\min}\\limits_{x \\in \\mathbb{R}^n}f(x), \\ \\mathrm{s.t.}\\ x \\in \\mathcal{X}\\ \\subset \\mathbb{R}^n \ ag{1}$ 与无约束问题不同，约束优化问题中自变量不能任意取值，这导致许多无约束优化算法不能直接使用。例如梯度法中沿着负梯度方向下降所得的点未必是可行点，要寻找的最优解处目标函数的梯度也不是零向量。这使得约束优化问题比无约束优化问题要复杂许多。放张很漂亮的等高线图如下：

考虑一种简单的情况，假设问题约束中仅含等式约束，即考虑问题： $\\mathop{\\min}\\limits_{x \\in \\mathbb{R}^n}f(x), \\ \\mathrm{s.t.}\\ c_i(x)=0 , i \\in \\mathcal{\\xi},x \\in \\mathbb{R}^n\ ag{2}$ 其中 $\\xi$ 为等式约束的指标集， $c_i(x)$ 为连续函数。

罚函数法的思想是将约束优化问题式(2) 转化为无约束优化问题来进行求解。为了保证解的逼近质量，无约束优化问题的目标函数为原约束优化问题的目标函数加上与约束函数有关的惩罚项。对于可行域外的点，惩罚项为正，即对该点进行惩罚；对于可行域内的点，惩罚项为0，即不做任何惩罚。因此，惩罚项会促使无约束优化问题的解落在可行域内。

等式约束的二次罚函数：对优化问题式(2)定义二次惩罚函数如下 $P_E(x,\\sigma)=f(x)+\\frac{1}{2}\\sigma \\sum\\limits_{i \\in \\xi}c_i(x),\\sigma > 0 \ ag{3}$ 这种罚函数对不满足约束的点进行惩罚，在迭代过程中点列一般处于可行域之外，因此它也被称为外点罚函数。

二次罚函数的特点如下：对于非可行点而言，当 $\\sigma$ 变大时，惩罚项在罚函数中的权重加大，对罚函数求极小，相当于迫使其极小点向可行域靠近；在可行域中，惩罚函数的全局极小点与约束最优化问题式(2) 的最优解相同。

举例：考虑如下优化问题及其对应的二次惩罚函数如下 $\\min \\ x+\\sqrt{3}y, \\ \\mathrm{s.t.}\\ x^2 + y ^2=1 \\Rightarrow P_E(x,y,\\sigma)=x+\\sqrt{3}y+\\frac{\\sigma}{2}(x^2+y^2-1)^2\ ag{4}$ 容易求得最优解为 $(-\\frac{1}{2},-\\frac{\\sqrt{3}}{2})$ ，下图给出了不同惩罚因子对用的二次罚函数的等高线图。

可以看出随惩罚因子的增大，罚函数的最小值与原问题的最小值越来越接近。但最优点附近的等高线越来越趋于扁平，这导致求解无约束优化问题的难度变大。在实际情况可能由于惩罚因子过小时，罚函数可能无下界，主要原因不可行点处的函数下降抵消了罚函数对约束违反的惩罚。实际上所有外点罚函数法均存在这个问题，因此惩罚因子的初值选取不应该太小。

二次罚函数算法的流程如下图所示。

不等式约束优化问题有如下形式： $\\mathop{\\min}\\limits_{x \\in \\mathbb{R}^n}f(x), \\ \\mathrm{s.t.}\\ c_i(x) \\leq 0 , i \\in \\mathcal{I},x \\in \\mathbb{R}^n\ ag{5}$ 不等式约束的二次罚函数: 只对 $c_i(x) >0$ 的点进行惩罚，定义如下： $P_I(x,\\sigma)=f(x)+\\frac{1}{2}\\sigma\\sum\\limits_{i \\in \\mathcal{I}}\ ilde{c}_i^2(x),\ ilde{c}_i(x)=\\max\\{c_i(x),0\\},\\sigma>0 \ ag{6}$ 注意上式可以利用梯度(可导)进行求解，但一般不是二阶可导的，不能直接利用二阶算法(牛顿法)求解子问题。

一般约束的二次罚函数：对于同时包含不等式与等式约束的优化问题，定义罚函数如下 $P_I(x,\\sigma)=f(x)+\\frac{1}{2}\\sigma[\\sum\\limits_{i \\in \\xi}c_i(x)+\\sum\\limits_{i \\in \\mathcal{I}}\ ilde{c}_i^2(x)],\\sigma>0 \ ag{7}$ 只需要将两种约束的罚函数相加就能得到一般约束优化问题的二次罚函数。

前面介绍的二次罚函数均属于外点罚函数，即在求解过程中允许自变量位于原问题可行域之外，当罚因子趋于无穷时，子问题最优解序列从可行域外部逼近最优解。自然地，如果想要使得子问题最优解序列从可行域内部逼近最优解，则需要构造内点罚函数。顾名思义，内点罚函数在迭代时始终要求自变量不能违反约束，因此它主要用于不等式约束优化问题。

考虑含不等式约束的优化问题式(5)，为了使得迭代点始终在可行域内，当迭代点趋于可行域边界时，我们需要罚函数趋于正无穷，常用的罚函数是对数罚函数。

对数罚函数：对于不等式约束最优化问题，定义对数罚函数如下： $P_I(x,\\sigma)=f(x)-\\sigma \\sum\\limits_{i \\in \\mathcal{I}}\\ln(-c_i(x)),\\mathbf{dom}P_I=\\{x|c_i(x)<0\\},\\sigma>0 \ ag{8}$ 对数罚函数的极小值严格位于可行域内部。原问题最优解通常位于可行域边界，此时需要调整惩罚因子使其趋近于0，这会减弱对数罚函数在边界附近的惩罚效果。

实例： 考虑如下优化问题及其对应的对数罚函数如下 $\\min x^2 +2xy+y^2+2x-2y,\\mathrm{s.t.}\\ x\\geq0,y\\geq 0 \\Rightarrow P_I(x,y,\\sigma)=x^2 +2xy+y^2+2x-2y -\\sigma(\\ln x + \\ln y) \ ag{9}$ 不同惩罚因子对应的罚函数的等高线如下图所示，可以看出，随着惩罚因子的减小，对数罚函数的最小值点和原问题最小值点越来越接近，但当变量趋于可行域边界时，对数罚函数趋于无穷大。

对数罚函数方法对应计算流程图如下所示。

和二次罚函数法不同，对数罚函数算法要求初始点为一个可行解，这是根据对数罚函数法本身的要求。常见的收敛准则可以为： $|\\sigma_k \\sum\\limits_{i \\in \\mathcal{I}}\\ln(-c_i(x^{k+1}))|\\leq \\epsilon,\\epsilon >0 \ ag{10}$ 实际以上算法的迭代点序列满足下式： $\\lim\\limits_{k \ o \\infty}\\sigma_k \\sum\\limits_{i \\in \\mathcal{I}}\\ln(-c_i(x^{k+1}))=0\ ag{11}$ 同样地，内点罚函数法也会有类似外点罚函数法的数值困难,即当惩罚因子趋于0时，罚函数的海森矩阵条件数趋于无穷大，因此子问题的求解将会越来越困难。

前面介绍的罚函数方法一个共同特点就是在求解的时候必须令罚因子趋于正无穷或零，这会带来一定的数值困难。而对于有些罚函数，在问题求解时不需要令罚因子趋于正无穷（或零），这种罚函数称为精确罚函数。

精确罚函数：罚因子选取适当，对罚函数进行极小化得到的解恰好就是原问题的精确解。这个性质在设计算法时非常有用，使用精确罚函数的算法通常会有比较好的性质。

?1罚函数：对应一般约束最优化问题，定义为 $P(x,\\sigma)=f(x)+\\sigma[\\sum\\limits_{i \\in \\xi}|c_i(x)|+\\sum\\limits_{i \\in \\mathcal{I}}\ ilde{c}_i^2(x)],\\sigma>0 \ ag{12}$ 对于精确罚函数，当罚因子充分大（不需要是正无穷）时，原问题的极小值点就是?1 罚函数的极小值点。

增广拉格朗日函数法的每一步构造一个增广拉格朗日函数，而该函数的构造依赖于拉格朗日函数和约束的二次罚函数。对于等式约束优化问题，增广拉格朗日函数定义为 $L_{\\sigma}(x,\\lambda)=f(x)+\\sum\\limits_{i \\in \\xi}\\lambda_ic_i(x)+\\frac{1}{2}\\sigma\\sum\\limits_{i \\in \\mathcal{\\xi}}\ ilde{c}_i^2(x),\\sigma>0 \ ag{13}$ 在拉格朗日函数的基础上，添加约束的二次罚函数。对于等式约束优化问题，最优解与最优乘子具有以下关系： $\ riangledown f(x^{\\ast})+\\sum\\limits_{i \\in \\xi}\\lambda_i^{\\ast}\ riangledown c_i(x^{\\ast})=0\ ag{14}$ 增广拉格朗日函数法计算流程如下图所示。

实例：考虑式(4)优化问题，对应的增广拉格朗日函数如下： $L_{\\sigma}(x,y,\\lambda)=x+\\sqrt{3}y+\\lambda(x^2+y^2-1)+\\frac{\\sigma}{2}(x^2+y^2-1)^2\ ag{15}$ 增广朗日函数法与二次罚函数法求解的点与最优点的位置关系如下图所示，可以得出前者得到解更接近于最优解，即说明增广朗日函数法相较于二次罚函数法在控制约束违反度上的优越性。

考虑一般约束优化问题如下： $\\min \\ f(x), \\mathrm{s.t.}\\ c_i(x)=0,i \\in \\mathcal{\\xi}, c_i(x)\\leq 0 ,i \\in \\mathcal{I}\ ag{16}$ 对于带不等式约束的优化问题，先通过引入松弛变量将不等式约束转化为等式约束和简单的非负约束，再对保留非负约束形式的拉格朗日函数添加等式约束的二次罚函数来构造增广拉格朗日函数。

通过引入松弛变量可以得到如下等价形式: $\\min \\ f(x), \\mathrm{s.t.}\\ c_i(x)=0,i \\in \\mathcal{\\xi}, c_i(x)+s_i=0 ,i \\in \\mathcal{I}, s_i \\geq 0 ,i \\in \\mathcal{I}\ ag{17}$ 保留非负约束，可以构造拉格朗日函数: $L(x,s,\\lambda,\\mu)=f(x)+\\sum\\limits_{i \\in \\xi}\\lambda_ic_i(x)+\\sum\\limits_{i \\in \\mathcal{I}}\\mu_i(c_i(x)+s_i) \ ag{18}$ 式(17)中等式约束的二次罚函数为 $p(x,s)=\\sum\\limits_{i \\in \\xi}c_i^2(x)+\\sum\\limits_{i \\in \\mathcal{I}}(c_i(x)+s_i)^2\ ag{19}$ 构造增广拉格朗日函数如下: $L_{\\sigma}(x,s,\\lambda,\\mu)=f(x)+\\sum\\limits_{i \\in \\xi}\\lambda_ic_i(x)+\\sum\\limits_{i \\in \\mathcal{I}}\\mu_i(c_i(x)+s_i)+\\frac{\\sigma}{2}p(x,s),s_i \\geq 0, i\\in \\mathcal{I},\\sigma \\geq 0 \ ag{20}$ 增广拉格朗日函数法: 通过增广拉格朗日函数，求解如下问题 $\\mathop{\\min}\\limits_{x,s}\\ L_{\\sigma_k}(x,s,\\lambda^k,\\mu^k),\\mathrm{s.t.}\\ s\\geq 0 \ ag{21}$ 求解上述问题的一个有效的方法是投影梯度法;另外一种方法是消去 $s$ ，求解只关于 $x$ 的优化问题，固定 $x$ ，关于 $s$ 的子问题可以表示为 $\\mathop{\\min}\\limits_{s}\\sum\\limits_{i \\in \\mathcal{I}}\\mu_i(c_i(x)+s_i)+\\frac{\\sigma_k}{2}\\sum\\limits_{i \\in \\mathcal{I}}(c_i(x)+s_i)^2\ ag{22}$ 根据凸优化问题的最优性理论， $s$ 为以上问题的一个全局最优解，当且仅当 $s_i=\\max \\{-\\frac{\\mu_i}{\\sigma_k}-c_i(x),0\\}, i \\in \\mathcal{I}\ ag{23}$ 将上式代入增广拉格朗日函数为 $L_{\\sigma_k}(x,\\lambda^k,\\mu^k)=f(x)+\\sum\\limits_{i \\in \\xi}\\lambda_ic_i(x)+\\frac{\\sigma_k}{2}\\sum\\limits_{i \\in \\xi}c_i^2(x)+\\frac{\\sigma_k}{2}\\sum\\limits_{i \\in \\mathcal{I}}(\\max \\{-\\frac{\\mu_i}{\\sigma_k}-c_i(x),0\\}^2-\\frac{\\mu_i^2}{\\sigma_k^2})\ ag{24}$ 其为关于x 的连续可微函数( $f(x),c_i(x),i \\in \\mathcal{I}\\cup \\xi$ 连续可微)，式(22)等价于 $\\mathop{\\min}\\limits_{x \\in \\mathbb{R}^n}\\ L_{\\sigma_k}(x,\\lambda^k,\\mu^k)\ ag{25}$ 并可以利用梯度法进行求解。

对于式(17)，其最优解需满足KKT条件如下： $0=\ riangledown f(x^\\ast)+\\sum\\limits_{i \\in \\xi}\\lambda_i^\\ast \ riangledown c_i(x^\\ast)+\\sum\\limits_{i \\in \\mathcal{I}}\\mu_i^\\ast \ riangledown c_i(x^\\ast), \\lambda_i^\\ast, \\mu_i^\\ast \\geq 0, i \\in \\mathcal{I}\ ag{26}$ 式(21)最优解 $x^{k+1},s^{k+1}$ 满足下式： $0=\ riangledown f(x^{k+1})+\\sum\\limits_{i \\in \\xi}(\\lambda_i^k+\\sigma_kc_i(x^{k+1}))\ riangledown c_i(x^{k+1})+\\sum\\limits_{i \\in \\mathcal{I}}(\\mu_i^k+\\sigma_k(c_i(x^{k+1})+s_i^{k+1}))\ riangledown c_i(x^{k+1})\ ag{27}$ $s_i^{k+1}=\\max\\{-\\frac{\\mu_k^k}{\\sigma_i^k},0\\}, i \\in \\mathcal{I}\ ag{28}$

可以计算乘子的更新公式如下： $\\lambda_i^{k+1}=\\lambda_i^k+\\sigma_kc_i(x^{k+1}), i \\in \\xi \ ag{29}$ $u_i^{k+1}=\\max \\{\\mu_i^k+\\sigma_k c_i(x^{k+1}),0 \\}, i \\in \\mathcal{I}\ ag{30}$

至此，给出约束优化问题的增广拉格朗日函数法及其对应的参数更新方式。

考虑凸优化问题如下： $\\mathop{\\min}\\limits_{x \\in \\mathbb{R}^n}\\ f(x), \\mathrm{s.t.}\\ c_i(x) \\leq 0, i=1,2,\\cdots,m \ ag{31}$ 其中 $f: \\mathbb{R}^n \\rightarrow \\mathbb{R}$ ， $c_i:\\mathbb{R}^n \\rightarrow \\mathbb{R}, i=1,2,\\cdots,m$ 为闭凸函数。

其增广拉格朗日函数为： $L_{\\sigma}(x,\\lambda)=f(x)+\\frac{\\sigma}{2}\\sum\\limits_{i=1}^{m}(\\max\\{\\frac{\\lambda_i}{\\sigma}+c_i(x),0\\}^2 -\\frac{\\lambda_i^2}{\\sigma^2})\ ag{32}$ 上式对应增广拉格朗日函数法更新公式为

$\\begin{cases}x^{k+1}\\approx \\mathop{\\mathrm{argmax}}\\limits_{x \\in \\mathbb{R}^n}\\ L_{\\sigma_k}(x,\\lambda^k) \\\\ \\lambda^{k+1}=\\lambda^k+\\sigma_k \ riangledown_{\\lambda}L_{\\sigma_k}(x^{k+1},\\lambda^k)=\\max\\{0, \\lambda^k+\\sigma_kc(x^{k+1}) \\}\\\\ \\sigma_k \鯋ow \\sigma_{\\infty}\\end{cases}\ ag{33}$ 式(31)的增广拉格朗日函数法计算流程如下所示。

考虑BP问题如下： $\\mathop{\\min}\\limits_{x \\in \\mathbb{R}^n}\\ ||x||_1, \\mathrm{s.t.}\\ Ax=b, A \\in \\mathbb{R}^{m\ imes n}(m \\leq n), b \\in \\mathbb{R}^m,x\\in \\mathbb{R}^n\ ag{34}$ 引入拉格朗日乘子，得到拉格朗日函数为 $L(x,y)=||x||_1+y^T(Ax-b), y\\in \\mathbb{R}^m\ ag{35}$ 其对偶函数为 $g(y)=\\mathop{\\inf}\\limits_x \\ L(x,y)=\\begin{cases}-b^Ty, ||A^Ty||_{\\infty}\\leq 1 \\\\ -\\infty, \\mathrm{others}\\end{cases}\ ag{36}$ 因此，得到如下对偶问题 $\\mathop{\\min}\\limits_{y \\in \\mathbb{R}^m}\\ b^Ty, \\mathrm{s.t.}||A^Ty||_{\\infty}\\leq 1 \ ag{37}$ 原始问题的增广拉格朗日函数法迭代更新公式如下：

$\\begin{cases}x^{k+1}=\\mathop{\\mathrm{argmin}}\\limits_{x \\in \\mathbb{R}^n}\\{||x||_1 +\\frac{\\sigma}{2}||Ax-b+\\frac{\\lambda^k}{\\sigma}||_2^2 \\}\\\\ \\lambda^{k+1}=\\lambda^k + \\sigma(Ax^{k+1}-b)\ ag{38}\\end{cases}$ 对于BP 问题，固定的σ 也可以保证增广拉格朗日函数法收敛,如下图所示。

引入松弛变量，原始问题的对偶问题式(37)可以等价写成： $\\mathop{\\min}\\limits_{y \\in \\mathbb{R}^m,s\\in \\mathbb{R}^n}\\ b^Ty, \\mathrm{s.t.}A^Ty-s=0,||s||_{\\infty}\\leq 1 \ ag{39}$ 引入拉格朗日乘子和罚因子，增广拉格朗日函数为 $L_{\\sigma}(y,s,\\lambda)=b^Ty+\\lambda^T(A^Ty-s)+\\frac{\\sigma}{2}||A^T-s||_2^2,||s||_{\\infty}\\leq 1\ ag{40}$ 其迭代公式如下所示： $\\begin{cases}(y^{k+1},s^{k+1})=\\mathop{\\mathrm{argmin}}\\limits_{y,||s||_{\\infty}\\leq1}\\{b^Ty+\\frac{\\sigma_k}{2}||A^Ty-s+\\frac{\\lambda^k}{\\sigma_k}||_2^2 \\}\\\\ \\lambda^k=\\lambda^k+\\sigma_k(A^Ty^{k+1}-s^{k+1})\\\\ \\sigma_{k+1}=\\min (\\rho \\sigma_k, \\bar \\sigma), \\rho >1, \\bar \\sigma < +\\infty \\end{cases}\ ag{41}$ 消去 $s$ 的增广增广拉格朗日函数法为:

$\\begin{cases}y^{k+1}=\\mathop{\\mathrm{argmin}}\\limits_{y}\\{ b^Ty +\\frac{\\sigma}{2}||\\Psi(A^Ty+\\frac{\\lambda^k}{\\sigma_k})||_2^2 \\}\\\\ \\lambda^k=\\sigma_k \\Psi(A^Ty^{k+1}+\\frac{\\lambda^k}{\\sigma_k})\\\\ \\sigma_{k+1}=\\min (\\rho \\sigma_k, \\bar \\sigma), \\rho >1, \\bar \\sigma < +\\infty \\\\ \\Psi(x)=\\mathrm{sign}(x) \\max \\{|x|-1,0 \\}\\\\ \\end{cases}\ ag{42}$ 其中： $\\mathrm{sign}(x)=\\begin{cases}1, x>0 \\\\ 0, x=0 \\\\ -1, x<0 \\end{cases}\ ag{43}$

半正定规划（Semidefinite programming，SDP）是凸优化问题的一个分支，它具有线性目标函数（由用户指定的最大化或最小化函数），且其定义在半正定矩阵构成的凸锥与仿射空间的交集上，即光谱面。

半正定规划问题如下： $\\mathop{\\min}\\limits_{X \\in S^n}<C,X>,\\mathrm{s.t.}<A_i,X≥b_i, i=1,2,\\cdots,m,X\\geq 0 \ ag{44}$ 其对偶问题为 $\\mathop{\\min}\\limits_{y \\in \\mathbb{R}^m}\\ -b^Ty, \\ \\mathrm{s.t.}\\sum\\limits_{i=1}^m y_iA_i \\leq C \ ag{45}$ 对于原始问题,增广拉格朗日函数为 $L_{\\sigma}=<C,X> -\\lambda^T(\\mathcal{A}(X)-b)+\\frac{\\sigma}{2}||\\mathcal{A}(X)-b||_2^2,X\\geq 0 \ ag{46}$ 其中： $\\mathcal{A}(X)=(<A_1,X>,<A_2,X>,\\cdots,<A_m,X>)^T\ ag{47}$ 那么，增广拉格朗日函数法为

$\\begin{cases}X^{k+1}\\approx \\mathop{\\mathrm{argmin}}\\limits_{X \\in S_+^n}L_{\\sigma_k}(X,\\lambda^k)\\\\ \\lambda^{k+1}=\\lambda^k-\\sigma_k(\\mathcal{A(X^{k+1})-b})\\\\ \\sigma_{k+1}=\\min (\\rho \\sigma_k, \\bar \\sigma), \\rho >1, \\bar \\sigma < +\\infty \\end{cases}\ ag{48}$ 其对偶问题对应的增广拉格朗日函数法为：

$\\begin{cases}y^{k+1}\\approx \\mathop{\\mathrm{arg \\min}}L_{\\sigma_k}(y,\\Gamma^k)\\\\ \\Gamma^{k+1}=\\sigma_k \\mathcal{P_{S_+^n}}(\\sum\\limits_{i=1}^my_i^{k+1}A_i-C+\\frac{\\Gamma^k}{\\sigma_k})\\\\ \\sigma_{k+1}=\\min (\\rho \\sigma_k, \\bar \\sigma), \\rho >1, \\bar \\sigma < +\\infty \\end{cases}\ ag{49}$ 其中 $\\mathcal{P_{S_+^n}}$ 为到半定锥集 $S_+^n$ 投影算子。

线性规划是非常经典的约束优化问题，它的目标函数和约束都是线性函数。因为其形式简单，在现实中有非常多的应用，线性规划一直受到人们的格外关注。求解线性规划问题的算法非常之多，最经典的要数Dantzig 在1947 年提出的单纯形法。由于线性规划问题具有特殊结构，它的解必然是在可行域的顶点（或某一边界处）取到，而单纯形法则是通过某种方式不断列出可行域的顶点然后一步一步寻找问题的最优解。由于线性规划可行域的顶点数可能多达 $\\mathbf{o}(2^n)$ 个（n 为自变量维数），因此单纯形法最坏情况下的复杂度是指数量级。实际上我们也可以构造出特殊的例子，使得单纯形法遍历可行域中的每一个顶点。这一现象表明对于某些大型问题和病态问题，单纯形法的效果可能很差，我们必须寻找其他办法来求解线性规划问题。

在大约30 年后，内点法应运而生，其中比较实用的算法是Karmarkar在1984 年提出的线性规划算法。内点法是在可行域内部寻找一条路径最终抵达其边界，这和单纯形法有着截然不同的思想。由于迭代点处于可行域内部，因此求解每个子问题的计算代价都远高于仅仅在可行域边界移动的单纯形法。然而内点法的一步迭代对问题解的改善是显著的，正因为如此，可以证明内点法实际上是一个多项式时间算法。

线性规划的原始问题和对偶问题分别如下所示： $(P) \\min c^Tx ,\\mathrm{s.t.}\\ Ax=b, x\\geq 0 \ ag{50}$