SCA凸优化实战：从非凸问题到迭代求解的完整指南

简介围绕SCA顺序凸逼近算法提供MATLAB平台下的凸优化实现代码适合正在学习凸优化理论、研究非凸问题求解以及从事信号处理、无线通信或能源系统优化等领域的工程师和研究人员阅读参考。SCA通过连续凸近似把非凸问题拆解为一系列易解的凸子问题是工程优化中非常实用的方法。压缩包共2个文件均为m脚本体积仅3KB其中sca.m是通用算法框架xiao_power_beizeng100.m则针对特定功率分配或通信场景给出完整示例便于对照理解算法结构与参数设置。目前已有2698人学习下载。读者可借助代码走通SCA的近似、求解、更新三个核心步骤学习如何用Taylor展开或松弛等方式构造凸近似函数并掌握MATLAB优化函数的调用与调试思路同时可根据示例快速迁移到自身任务减少重复开发成本。1. SCA凸优化在解决什么非凸问题也能迭代求解SCA凸优化不是某个现成工具箱的名字而是把非凸优化问题拆成一串凸子问题逐个求解的工程范式。只要目标函数或约束里有一个非凸项CVXPY这类凸优化求解器就会直接报 DCP error而 SCA 的思路是在当前点附近把非凸部分替换成凸的代理函数解出一个子问题拿到新点再重复直到收敛。这套做法在无线功率分配、传感器定位、波束成形、稀疏恢复里几乎是标配。这篇文章从 SCA 的数学骨架写到可复现的 Python 算例把参数、停止条件、避坑和验证方法一次讲透适合已经具备基础优化知识、想拿 SCA 解决自己非凸问题的工程师和研究生。2. SCA的数学框架三种凸近似做法与收敛性判断2.1 为什么非凸问题直接让CVXPY哑火凸优化求解器能工作的前提是问题满足 DCPDisciplined Convex Programming规则目标函数和约束里的每一项都得能被识别为凸或仿射。工程里的大多数非凸问题问题本身并不大难在表达式结构上。一个典型的例子是测距定位的极大似然估计[ \min_{x} \sum_{i1}^{N} (|x - a_i| - d_i)^2 ]其中 (a_i) 是已知基站坐标(d_i) 是带噪测距。单独看 (|x - a_i|) 是凸函数减一个常数再平方听起来好像还是凸的。但注意平方这个外层函数 ((t - d_i)^2) 在 (t) 较大时单调上升、较小时单调下降复合后不再满足凸函数复合规则。整项不凸CVXPY 会直接告诉你 “Problem does not follow DCP rules”。这类问题还有个共同点非凸性常常只来自少数几项其余全是凸的。如果能把非凸项在当前迭代点附近换成一个凸的近似表达式让整个子问题变成凸问题就可以继续用 CVXPY、OSQP、MOSEK 这类成熟求解器。SCA 的核心就在这个“换”字上怎么换、换完怎么保证迭代不跑飞、收敛到哪里是三个真正要解决的问题。2.2 SCA的迭代骨架与三种凸化方式SCA 的迭代骨架非常统一。假设原始问题是 (\min f(x) \text{ s.t. } x \in \mathcal{C})其中 (f) 非凸、(\mathcal{C}) 是凸可行域。在第 (k) 步算法在当前点 (x_k) 附近构造一个凸代理函数 (\tilde{f}(x; x_k))求解凸子问题[ x_{k1} \arg\min_{x \in \mathcal{C}} \tilde{f}(x; x_k) ]然后更新 (x_k)重复直到收敛。工程上用到的主要有三种凸化方式各有代价。第一种是一阶泰勒展开也是使用频率最高的。把非凸光滑项在当前点展开成线性项因为线性函数天然仿射保留其余凸结构不变就可以凑出一个凸子问题。代价是近似只在一个小邻域内有效所以通常会配合信任域约束或阻尼步长。第二种是二次上界也叫 majorization。若 (f) 的梯度满足 Lipschitz 连续常数为 (L)可以构造[ \tilde{f}(x; x_k) f(x_k) \nabla f(x_k)^T (x - x_k) \frac{L}{2} |x - x_k|^2 ]这个代理函数是凸的并且全局位于原函数上方。它的好处是理论上可以保证目标值单调下降坏处是 (L) 估计不准时步长和收敛速度都会受影响。第三种是变量代换。举个例子通信波束成形里经常出现 (|w^H h|^2) 这类项通过引入辅助变量和相位旋转把非凸的二次型转成一组凸约束。这不是 SCA 专属但经常和 SCA 搭配使用先用代换把问题“预处理”得尽量凸剩下少量顽固非凸项再用前两种方法处理。选择哪种方式主要看问题结构。一个实用的判断标准是如果非凸项是光滑的优先一阶泰勒加信任域如果不光滑但梯度有界用二次上界如果能通过代换把非凸结构消掉先消掉再做 SCA近似出来的子问题更紧。2.3 收敛性判断上界逼近与单调下降SCA 的收敛性在数学上可以走 MMMajorization-Minimization算法的框架。只要代理函数满足三个条件在 (x_k) 处与原函数同值、一阶条件相同、是原函数的全局上界或至少在信任域内成立那么每次迭代后目标值不会上升。这个单调性保证了迭代序列不会发散多个聚点满足原问题的 KKT 条件。但这里要泼一盆冷水SCA 提供的是“稳定点”不是“全局最优点”。对非凸问题SCA 和梯度下降一样会卡进局部最优。工程上真正的收敛性判断要回答两个问题一是迭代有没有真的收敛到不动点二是这个不动点是不是够好。第一个看停止准则第二个看多初值对照这两个话题分别在第四章和第六章展开。一个常见误区是把“目标值下降”当成“找到最优解”。目标值单调下降只能说明代理函数构造得不错不代表解的质量高。很多论文里的收敛性证明也只覆盖到“收敛到稳定点”工程落地时必须自己对解做验证否则结果可能看着漂亮放到真实系统里误差很大。下表是三种凸化方式的直接对比便于在拿到新问题时快速选型凸化方式适用对象代理函数特征工程代价一阶泰勒 信任域光滑非凸项、变量维度高线性近似局部有效需要调信任域半径容易振荡二次上界majorization梯度 Lipschitz 连续的目标全局上界收敛稳需要估计 Lipschitz 常数偏保守变量代换预凸化二次型、绝对值结构表达式改写不损失精度依赖问题结构通用性有限3. 用PythonCVXPY跑通SCA传感器定位的完整算例3.1 问题设定为什么测距定位的目标函数是非凸的这里用第二章的测距定位例子完整走一遍。假设平面上有 5 个基站坐标已知目标节点在某个位置 (x)我们能拿到带噪声的测距 (d_i)。最直接的做法是让估计位置对应的距离尽量贴近实测距离这个最小二乘问题虽然在工程上用了很多年但它本质上不是凸问题尤其在测距噪声较大时直接用凸优化求解器是跑不过去的。SCA 的做法是把目标函数里的 (|x - a_i|) 在当前位置 (x_k) 处做一阶泰勒展开。注意这里的技巧(|x - a_i|) 对 (x) 的梯度是 (\frac{x - a_i}{|x - a_i|})是沿基站指向目标方向的单位向量。把它展开成线性项后残差部分变成关于 (x) 的仿射函数平方后就是凸二次函数整个子问题就变成了一个带信任域约束的二次规划。基线选择上可以用基站坐标平均值作为初值。这个初值不一定准但能让第一轮迭代不至于偏离太远之后再靠 SCA 迭代修正。3.2 代码实现一次SCA迭代中的凸化与求解下面的代码用 NumPy 生成仿真数据用 CVXPY 在每个 SCA 迭代中求解凸子问题。核心是每次迭代都重新构造一次代理函数并用信任域约束限制 (x) 的移动范围。import numpy as np import cvxpy as cp # ---- 生成仿真数据 ---- np.random.seed(7) A np.array([ [0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0], [5.0, 13.0] ]) # 5 个基站坐标 x_true np.array([4.2, 3.8]) # 目标真实位置 dist_true np.linalg.norm(A - x_true, axis1) d dist_true np.random.normal(0, 0.05, A.shape[0]) # ---- SCA 参数 ---- x_cur np.array([5.0, 5.0]) # 初始点基站坐标均值 sigma 1.0 # 信任域半径 max_iter 30 eps 1e-5 # 停止阈值 obj_hist [] # ---- SCA 主循环 ---- for it in range(max_iter): # 当前点到各基站的距离及单位方向向量 diff x_cur - A # (N, 2) dist np.linalg.norm(diff, axis1) grad diff / dist[:, None] # 距离函数在当前点的梯度方向 # 定义优化变量 x cp.Variable(2) # 一阶泰勒展开代替 ||x - a_i|| linear_dist dist grad (x - x_cur) # 凸近似目标残差平方和 residual linear_dist - d cost cp.sum_squares(residual) # 信任域约束 cons [cp.norm(x - x_cur) sigma] prob cp.Problem(cp.Minimize(cost), cons) prob.solve() x_next x.value obj_hist.append(prob.value) # 停止条件位置变化足够小 if np.linalg.norm(x_next - x_cur) eps: break x_cur x_next print(SCA 估计位置:, x_cur) print(真实位置: , x_true) print(定位误差: , np.linalg.norm(x_cur - x_true))这段代码里最关键的是linear_dist dist grad (x - x_cur)。dist和grad都是 NumPy 常量数组x - x_cur是 CVXPY 变量表达式两者做矩阵乘法后得到一个 N 维仿射表达式。这个表达式就是对 (|x - a_i|) 在 (x_cur) 处的一阶近似残差平方和因此变成凸二次函数。三个参数直接影响迭代行为sigma是每次迭代允许的最大移动距离设置太大会让近似失效设置太小则收敛极慢eps控制停止精度max_iter是安全阀防止不收敛时无限循环。代码里没有直接显示目标函数值曲线建议自己把obj_hist画出来观察是不是单调下降这是判断 SCA 实现是否正确的最快方式。3.3 初始点与信任域从翻车到稳收敛的两个改动第一次跑上面代码时最常见的翻车点是初始点选得离真值太远。如果把x_cur从[5,5]改成[50,50]第一轮迭代的线性近似离真解很远信任域约束半径又限制每步只能走 1收敛步数会明显增加。更糟的情况是目标函数有多模态初值不对会把迭代拉进错误的谷里。工程上有两个补丁。第一个是初值不要拍脑袋先用线性最小二乘粗定位算出初始点。把 (|x - a_i|) 两边平方并展开可以整理成关于 (x) 的线性方程组用几行 NumPy 就能解出一个不错的起点。第二个是信任域自适应如果连续两轮目标值都在下降说明近似质量不错可以适当放大sigma一旦目标值上升立刻把sigma缩小一半。# 可选的线性粗定位初值 b 0.5 * (np.sum(A**2, axis1) - d**2) M np.hstack([A, np.ones((A.shape[0], 1))]) x_ls, _, _, _ np.linalg.lstsq(M, b, rcondNone) x_cur x_ls[:2]这段代码从测距方程平方展开出发构造线性方程组再用最小二乘解出初始位置。它的精度在高噪声下有限但作为 SCA 的起点远比均值初值可靠能让后续迭代更容易收敛到好解。4. 可调参数与停止条件让SCA收敛更快的五个旋钮4.1 代理函数的紧度一阶泰勒以外的选择一阶泰勒近似是所有凸化方式里最简单的一种但它不是免费的近似误差会随着 (|x - x_{k}|) 增大而变大。信任域约束相当于人为限制误差范围代价是每轮迭代只能走一小步。想走大步就要换更紧的代理函数。常用的一种替代是二阶近似在泰勒展开里加入 Hessian 项。二阶代理更贴近原函数收敛所需的迭代次数更少但 Hessian 不一定是半正定的需要做特征值修正或加正则项否则凸性无法保证。实际应用中我不会一上来就用二阶只有当一阶版本振荡太严重、收敛太慢时才升级。另一个思路是在代理目标里加入近端正则项比如把子问题改成[ \min_{x} \tilde{f}(x; x_k) \frac{\mu}{2} |x - x_k|^2 ]这个正则项的作用和信任域类似把 (x) 拉向当前点防止代理函数离当前点太远。和信任域相比它的优势是可以写成无约束形式方便套用已有的快速 QP 求解器劣势是 (\mu) 多了一个需要调节的参数。4.2 五个直接决定收敛快慢的参数SCA 的收敛速度不是算法决定的而是参数决定的。同一套代码参数不同可能是十步收敛也可能是一百步振荡。我把工程里最有影响的五个旋钮列出来参数典型范围调大效果调小效果注意事项信任域半径 ( \sigma )问题尺度的 0.1% ~ 10%迭代更快易振荡更稳收敛变慢位置问题看坐标尺度功率问题看对数域阻尼系数 ( \gamma )0.1 ~ 1.0步长更大可能越过好点更稳迭代增加用 (x_{k1}x_k\gamma(\hat{x}-x_k))正则化强度 ( \mu )0 ~ 10稳定但偏向旧点近似作用增强与信任域二选一即可初始点质量尽量靠近好解收敛快且结果稳可能进入错误局部解用线性近似或粗搜索给初值停止容差 ( \epsilon )(10^{-6} \sim 10^{-3})多算几轮结果更精提前停止结果偏粗只靠位移判据不够配合目标值变化信任域半径和阻尼系数是互相替代的关系。我一般会先固定一个只调另一个。如果两个都调参数空间变大出了问题很难定位是哪一项导致的振荡。sigma的初值可以按问题尺度设基站坐标跨度是 10sigma1.0起步就合理如果问题里的变量是功率取值范围是 0 到 1 瓦sigma就得用 0.01 这个量级。4.3 停止条件怎么设目标值、位移、梯度范数很多 SCA 实现只用位置变化量做停止条件也就是上面代码里的eps。这有一个盲区当 SCA 进入一个很平缓的区域位置每次只动一点点但离真正的稳定点还很远代码会提前退出。反过来目标值接近收敛时位置可能仍在小幅抖动只靠位移判断又会多跑很多轮。稳妥的停止条件至少有两个一是位置变化的二范数小于容差二是目标值相对变化小于容差两个条件同时满足才停。更严格的做法是检查一阶最优性残差也就是计算 (|x - \text{Proj}_{\mathcal{C}}(x - \nabla f(x))|)。这个值越小说明越接近 KKT 条件。对于没有约束的问题梯度范数就是现成判据对于带信任域约束的问题需要把投影算子算出来。实际工程里我不会把这套判据一开始就做全而是先跑通最小实现确认迭代正常下降后再加判据。过早追求完美停止条件反而会被参数调试拖住。5. SCA避坑清单五次实验里最常踩的五个问题5.1 目标值一直降但位置基本不动现象画出obj_hist目标函数单调下降看起来一切正常打印位置却发现迭代十轮几乎没挪动定位误差也没变小。原因这是信任域半径太小或阻尼系数过低的典型症状。代理函数被限制在一个很小的邻域内每轮只能往正确方向挪一点点目标值虽然下降但位置收敛极慢。另一种可能是梯度方向写反了——距离函数 (|x-a_i|) 的梯度是从基站指向目标点如果代码里写成了(A - x_cur) / dist方向正好相反迭代会在错误方向上反复试探。解决先检查梯度方向确认diff x_cur - A而不是A - x_cur。方向没问题就把sigma调大一倍或者把阻尼系数gamma提到 0.8 以上观察两三轮内位置移动量是否变大。5.2 凸化后的子问题仍然通不过DCP校验现象prob.solve()直接抛出异常提示问题不是 DCP。原因最常见的是把变量写进了非仿射表达式。比如想计算线性化距离时忘了dist和grad应该作为预先算好的常数把它们写成了依赖变量的表达式导致 CVXPY 无法识别凸性。另一个常见操作是尝试对变量计算sqrt(cp.sum_squares(x - A))这会让整项变成非凸。解决记住一个规则——SCA 的每一次迭代里原问题的非凸项必须完全被常数系数替代。代码里的grad和dist要在进入 CVXPY 之前用 NumPy 算好CVXPY 表达式中只允许出现grad (x - x_cur)这种仿射形式。如果还是要写cp.norm(x - a)说明这个非凸项没有被真正凸化。5.3 迭代中突然出现NaN或inf现象某些轮次计算出的x_next变成 NaN程序直接崩溃。原因几乎都是距离为 0 导致的除以 0。当目标节点恰好在某个基站附近时np.linalg.norm(diff, axis1)会算出 0grad diff / dist[:, None]产生无限值后续矩阵运算全部污染。这类问题在真实数据里尤其阴险因为仿真时不一定踩得到。解决算距离时加一个保护下限dist np.maximum(dist, 1e-8)再算梯度。这个操作不会影响正常数据的精度但在数值上避免了除零。如果问题里存在多个基站还需要检查有没有波峰波谷导致的极端值必要时对测距数据加截断。5.4 SCA收敛了但结果明显偏离真值现象目标值正常下降、位置也稳定在某个点把这个点代回去算定位误差发现离真值很远。原因SCA 收敛到了局部最优。非凸问题的目标函数常常有多个谷底SCA 作为一种局部方法没有能力判断自己落在哪个谷里。初始点选得不好就会稳定地收敛到错误的位置。解决多初值启动是最直接的后悔药。固定其它参数不变取 5 个不同初值分别跑 SCA比较每个结果的最终目标值。如果多个初值收敛到不同位置选目标值最低的那个如果多个初值都收敛到同一位置结果可信度就高很多。这个方法算下来只多几倍计算量但能避开大部分局部最优问题。5.5 把非凸等式约束照抄进子问题现象想在定位问题里加入 (|x-a_i| d_i) 的约束来提高精度结果 CVXPY 报错或者求出来的解直接不可行。原因等距约束是典型的非凸约束。在 SCA 框架里目标函数的非凸项可以凸化但约束里的非凸项不能简单保留凸子问题里出现非凸可行域会让整个求解过程失效。解决两条路。一是把等式约束改成不等式约束(|x-a_i| \le d_i \delta)这个约束在 (x) 上是凸的二是彻底放弃硬约束目标函数里把残差项权重调大用惩罚代替约束。工程上第二种方法更常用因为噪声数据下硬约束本身就不合理惩罚项可以直接用 (w_i (|x-a_i| - d_i)^2) 的形式权重 (w_i) 按测距噪声方差设置。6. 验证SCA求解质量的三个技巧对偶间隙、KKT检验与多初值很多人是从凸优化的对偶理论和 KKT 条件那几页 PPT 入门这门课的但到 SCA 落地时这些理论不是用来背定义的而是用来验证结果的。第一个技巧是看子问题的对偶间隙。每个 SCA 迭代都是一次凸优化求解器如果不报错原始问题和对偶问题的最优值间隙应该趋近于零。如果求解器的状态不是 optimal说明子问题本身的数值稳定性有问题这时应该回头查参数而不是继续迭代。第二个技巧是检查 KKT 残余。SCA 收敛后把最终点 (x^) 代回原始目标函数计算梯度投影残差 (|x^- \text{Proj}_{\mathcal{C}}(x^* - \nabla f(x^*))|)。这个值越小说明越接近一阶最优性条件。这个检查不依赖求解器内部的报告完全可以用 NumPy 自己算作为最终解质量的硬指标。第三个技巧是多初值一致性验证。把第三章的 SCA 循环封成run_sca(A, d, x0)函数后分别用几个差异明显的初值跑比较最终位置和最终目标值starts [np.array([2.0, 2.0]), np.array([8.0, 8.0]), np.array([5.0, 9.0])] for x0 in starts: x_opt run_sca(A, d, x0) print(初值, x0, - 估计, x_opt, 目标值, obj_value(x_opt))多个初值若收敛到同一目的位置这个解落在该局部区域的可能性很大若分散在不同位置选目标值最低的一个并标注该结果受初值影响需要进一步加约束或加正则化。我现在做 SCA 类方案的习惯是先跑最小实现再上参数自适应最后一定补多初值检查。这套流程在好几个项目里救过我最典型的一次是功率分配问题里单初值收敛很快但性能差了 3dB 以上换多初值后找到了明显更好的点。希望这些经验帮到你尤其是第一次在 CVXPY 里实现 SCA 时少走几步弯路。本文还有配套的精品资源点击获取