☰
SCA凸优化实战:非凸问题迭代逼近与MATLAB代码解析
2026/9/26 8:45:01 网站建设 项目流程

简介:一份专注SCA(Sequential Convex Approximation)凸优化算法的MATLAB实现资源包,面向需要处理非凸优化问题的学习者、研究者与工程技术人员,适用于无线通信、信号处理、能源系统等典型应用场景。压缩包内共2个文件,均为.m脚本,整体大小仅3KB,轻量精炼:一个提供SCA算法通用迭代框架,另一个给出特定问题的近似示例与调用方式,便于读者对照理解“凸近似—凸子问题求解—变量更新”的核心流程。内容还涉及凸函数与凸集判定、凸优化问题标准形式、Taylor展开近似构造技巧及收敛性分析要点,帮助读者从原理到代码掌握SCA方法。资源已有2698人学习下载,研读后可快速获得可直接运行的MATLAB脚本,并理解如何将非凸原问题分解为系列可解的凸子问题,避免从零搭建框架,适合入门学习与工程实践参考。

1. SCA 凸优化:非凸问题并非无解,关键是找对逼近路径

SCA 凸优化(Sequential Convex Approximation,顺序凸近似)在工程优化里出现的频率比想象得高,尤其在无线通信功率分配、波束成形、信号处理和能源系统调度这几类场景中。反直觉的一点是:很多人以为非凸问题只能靠随机搜索或启发式算法碰运气,但 SCA 的做法是把原问题按迭代拆成一串凸子问题,每一步用成熟的凸优化工具求解,最终收敛到一个满足 KKT 条件的驻点或满意解。这份 sca.zip 里只有两个文件——sca.m 和 xiao_power_beizeng100.m,一个负责通用迭代框架,一个负责演示具体功率分配问题的建模与求解。适合谁?如果你正在复现论文里的迭代算法,或者手头有一个能写出数学模型但调不通的非凸优化问题,这份代码能给你一条可以直接改着用的路径。

2. 先立住优化模型:凸函数判定、非凸来源与 SCA 三步迭代

2.1 凸函数与凸集:先搞清“局部最优等于全局最优”的成立条件

凸优化教科书里最常被引用的一句话是:凸问题上的任何局部最优解都是全局最优解。但这个结论成立有两个前提,一是目标函数是凸函数,二是可行集是凸集,两个条件缺一不可。凸函数的定义是:对任意两个点 x1、x2,以及任意 λ∈[0,1],都满足 f(λx1+(1−λ)x2) ≤ λf(x1)+(1−λ)f(x2)。这个不等式看起来抽象,落到图像上就是函数曲线上的任意两点连线,不会落在曲线下方。

工程上判断一个多维函数是否为凸函数,常用的做法是看它的 Hessian 矩阵是否半正定。如果 Hessian 矩阵的所有特征值都大于等于零,那函数就是凸的;如果特征值里有负数,函数就是非凸的。很多人只画一维曲线看“是不是碗形”,这在一维情形下没问题,但到了多维就不可靠了。我见过不少把非凸问题当成凸问题直接丢给 CVX 的案例,结果求解器要么报错,要么给出的结果明显不合理。

可行集的凸性同样重要。两个可行点之间的连线上的任意一点,也必须满足所有约束条件,这个集合才是凸集。线性约束自然满足这个性质,但非线性约束就不一定了,比如二次等式约束 h(x) = ‖x‖² = c 这样的集合就是一个球面,球面上两个点的连线会穿过球体内部,并不完全落在球面上,所以它是非凸集合。 SCA 的切入点就在这里:每轮迭代里用一个凸集合去近似这个非凸集合,让凸优化工具能正常工作。

2.2 非凸项的三个典型来源:分式、乘积与指数嵌套

实际工程问题里的非凸项并不是程序员故意制造出来的,它们往往来自物理模型的表达方式。第一个典型来源是分式项。通信领域里的 SINR(信干噪比)就是最典型的分式:分子是用户自身信号功率,分母是干扰加噪声,这个结构在优化变量上天然是非凸的。第二个典型来源是变量乘积,比如两个优化变量相乘 x·y,或者向量与自身共轭转置构成的外积 wwᴴ。功率分配问题里经常出现变量乘积后还要约束秩为一的情况,这种约束几乎都是非凸的。第三个典型来源是指数或对数嵌套,比如能效优化里的 log(1+SINR) 一旦和分母上的功率相加耦合,函数整体就失去了凸性。

这些非凸项不会因为我们“希望它是凸的”就自动变凸。处理它们有两个方向:一是做变量替换,把非凸结构转化成凸结构,但很多时候变量替换之后又会冒出新约束,等于按下葫芦浮起瓢;二是用 SCA,在每一轮迭代中把非凸项替换成一个在当前展开点附近成立的凸近似。SCA 的核心思想可以概括成一句话:不是直接解原问题,而是解一列逐渐逼近原问题的凸子问题。

2.3 SCA 的三步循环:近似、求解、更新

标准 SCA 每一步迭代都包含三个阶段。第一步是近似,在当前的迭代点 x(k) 处,把目标函数里非凸的部分替换成凸近似函数,同时对非凸约束做同样的处理。第二步是求解,用内点法、梯度投影法或者现成工具比如 MATLAB 的 quadprog、CVX 求解这个凸子问题,得到临时解 x_temp。第三步是更新,把 x_temp 作为下一轮迭代的展开点,也可以加入阻尼系数避免震荡。

整个循环的数学形式大致如下:

初始化 x(0) for k = 0, 1, 2, ...: 构造非凸项的凸近似 f_tilde(x; x(k)) 求解凸子问题: min f_tilde(x; x(k)) s.t. 凸约束集合 得到临时解 x_temp 更新: x(k+1) = x_temp (或带阻尼的更新) 判断是否收敛 end

这里每一步都有讲究。近似构造的质量直接决定收敛速度,如果近似函数偏离原函数太远,迭代会来回跳动甚至发散。求解阶段必须保证子问题严格凸,否则又落回局部最优的老问题。更新阶段如果直接让 x(k+1) = x_temp,在问题条件数很差时容易震荡,常见的处理是引入一个步长因子 α∈(0,1],把更新改成 x(k+1) = x(k) + α·(x_temp − x(k))。

提示:判断 SCA 是否走对方向,最直观的方式是看每次迭代后原目标函数值是否呈单调下降或单调上升趋势。虽然 SCA 并不保证每轮都严格单调,但在绝大多数设计良好的实现里,目标值曲线应该像一条台阶式下降的曲线,如果看到锯齿状震荡,优先怀疑步长和近似质量。

3. 拆解 MATLAB 代码:sca.m 的框架与 xiao_power_beizeng100.m 的功率分配示例

3.1 从文件名看资源结构:两个文件到底怎么分工

sca.zip 里只有两个 .m 文件,这种精简结构在论文复现包里很常见。sca.m 一看就是主程序或通用函数名,它承担的是 SCA 迭代框架:初始化、循环、停止判断、结果输出。xiao_power_beizeng100.m 从命名习惯看,xiao 是“小”,power 是“功率”,beizeng100 是“倍增100次”的意思,大概率是一个小规模的功率增长或功率分配演示脚本,循环 100 次,每次让某个功率量做倍增或者步进更新,用来观察 SCA 在动态场景下能不能跟踪最优解。

这两个文件的设计思路一般是:xiao_power_beizeng100.m 调用 sca.m,前者定义问题模型,后者执行 SCA 迭代。如果直接运行 xiao_power_beizeng100.m 就能出结果,说明它已经把问题相关的参数都内置了。如果你要换自己的问题,只需要保留 sca.m 的循环骨架,重写近似构造和求解那两段。

3.2 sca.m 通用框架:核心循环与停止条件

一个可复用的 sca.m 不会把具体的数学函数写死在代码里,而是把“近似构造”“求解”“目标值计算”分别写成子函数或函数句柄,这样换问题时不需要动主循环。常见的框架拆成三个函数块,具体写成 MATLAB 大概是这个样子:

function [x_opt, obj_hist] = sca(x0, approx_func, solve_func, obj_func, params) % sca.m -- SCA 通用迭代框架 % 输入: % x0 : 初始点(列向量) % approx_func : 函数句柄, 用于在给定点处构造凸近似模型 % solve_func : 函数句柄, 用于求解当前的凸子问题 % obj_func : 函数句柄, 用于计算原问题的真实目标值 % params : 结构体, 包含 max_iter, tol, alpha 等参数 % 输出: % x_opt : 收敛后的最优解 % obj_hist : 每次迭代的原目标函数值, 用于画收敛曲线 x = x0(:); n = length(x); obj_hist = zeros(params.max_iter, 1); for k = 1:params.max_iter % 1. 近似:在当前迭代点 x 处构造凸子问题 model = approx_func(x); % 2. 求解:把凸子问题交给求解器 x_temp = solve_func(model, x); % 3. 更新:加入阻尼系数, 防止迭代点来回跳 alpha = params.alpha; x = x + alpha * (x_temp - x); % 记录真实目标值, 注意这里必须用原目标函数计算 obj_hist(k) = obj_func(x); % 停止条件: 相邻两次迭代的变量变化足够小 if norm(x_temp - x, inf) < params.tol obj_hist = obj_hist(1:k); break; end end x_opt = x; if k == params.max_iter fprintf('[sca] 达到最大迭代次数 %d, 未完全收敛\n', params.max_iter); end end

这段代码的逻辑很直接:外层 for 循环就是 SCA 三步循环,approx_func 和 solve_func 是两个函数句柄,这是 MATLAB 里做模块解耦最常用的方式。alpha 这个参数值得细说,它通常取 0.5 到 1 之间。alpha = 1 时就是完全信任凸子问题的解,收敛快但容易震荡;alpha 偏小时每步走得保守,曲线平滑但迭代次数变多。tol 一般取 1e-6 或 1e-8,如果你的变量本身量级很大,相对判断会比绝对判断更稳。

这里还有一个容易被忽略的细节:obj_hist 记录的一定要用原问题的目标函数,而不是近似函数的目标值。近似函数每轮都在变,它的下降不代表原目标下降,只有原目标函数值才是判断收敛的可靠依据。我把这条写进过很多次代码注释里,因为确实有朋友拿近似目标值画收敛曲线,画出来一条漂亮的单调曲线,但实际问题根本没收敛。

3.3 xiao_power_beizeng100.m 的功率分配示例:变量、约束与监视量

xiao_power_beizeng100.m 如果按“小功率倍增 100 次”来理解,那它大概率是一个这样的演示:初始给每个用户一个很小的功率,循环 100 次让功率以固定倍数或步长增长,每一步调用 SCA 更新出一个可行功率向量,同时记录系统吞吐量或 SINR。这类脚本在无线通信论文复现里出现频率很高,它的核心结构往往长这样:

% xiao_power_beizeng100.m % 演示 SCA 在小规模功率分配问题上的迭代行为 % 场景: N 个用户共享一个信道, 每个用户有一个功率变量 p_i clear; clc; rng(1); N = 4; % 用户数 P_max = 1; % 归一化总功率上限 p0 = 0.01 * ones(N,1); % 每个用户初始小功率 % 构造信道增益矩阵, 对角线为用户自身信道, 非对角线为干扰信道 H = abs(randn(N) * 0.5 + 0.5); H(1:N+1:end) = 1.0; % 迭代参数 max_iter = 100; params.alpha = 0.8; params.tol = 1e-6; % 记录吞吐量和 SINR throughput_hist = zeros(max_iter, 1); sinr_hist = zeros(max_iter, N); p = p0; for k = 1:max_iter % 固定当前点, 计算每个用户的干扰功率 interference = H * p - diag(H) .* p; % 用当前点构造 SINR 约束的凸近似(分母线性化) sinr_threshold = 0.5; % 把 sinr >= threshold 改写成 信号 >= threshold * 干扰 % 这是一个线性约束, 是凸的 A_ineq = threshold * H; A_ineq(1:N+1:end) = A_ineq(1:N+1:end) - diag(H); b_ineq = zeros(N,1); % 用 quadprog 求解最小化 -sum(log(1+sinr)) 的凸子问题 % 这里简化为最小化 -sum(log(信号项)) H_diag = diag(H); f = -log(H_diag); % 线性目标系数 options = optimoptions('quadprog', 'Display', 'off'); p_new = quadprog(2*eye(N), f, A_ineq, b_ineq, ... ones(1,N), P_max, zeros(N,1), ones(N,1), p, options); % 阻尼更新 p = p + params.alpha * (p_new - p); % 计算真实 SINR 和吞吐量 sinr = (H_diag .* p) ./ (interference + 1e-6); throughput_hist(k) = sum(log2(1 + sinr)); sinr_hist(k,:) = sinr'; end

这段代码里的核心技巧是 SINR 约束的线性化:SINR 大于等于阈值这个非凸约束,在固定干扰项后可以改写成信号强度大于等于阈值乘干扰,这一步做完约束就变成了线性约束。quadprog 专门求解二次规划问题,这里因为目标函数被简化成了线性形式,所以代价矩阵用了 2*eye(N),实际项目中你要把真实的目标二阶信息填进去。A_ineq 矩阵的构造是整个例子的核心,对角线减去 diag(H) 是为了把自身信号从干扰项里扣除,这部分容易写错,建议写完之后打印出来逐个检查。

注意演示脚本里用了 rng(1) 固定随机种子,这意味着每次运行结果可复现。实际调参时我会盯两个监视量:throughput_hist 的曲线形态和 p 的逐维变化幅度。如果 throughput 曲线稳步上升最后趋于平缓,说明 SCA 的近似方向是对的;如果曲线突然跳降,大概率是某个约束线性化出了错,或者阻尼系数 alpha 太大导致迭代点越过了可行域边界。

4. 凸近似怎么构造:Taylor 展开、约束松弛与收敛判据

4.1 一阶 Taylor 展开:构造凸下界的最直接方式

SCA 里最常用的近似工具是一阶 Taylor 展开。对于目标函数中的非凸项 g(x),在展开点 x(k) 处做一阶展开,得到 g(x) ≈ g(x(k)) + ∇g(x(k))ᵀ(x − x(k)),这样就把原函数替换成仿射函数,而仿射函数既是凸函数也是凹函数,可以自由决定它作为上界还是下界。如果是凸函数被下放,一阶展开天然就是原函数的全局下界,这是凸函数的经典性质。

对于非凸项,一阶展开不一定是下界,也可能是上界,取决于具体函数结构。工程上常见的处理方法是:把目标函数分解成“凸部分 + 非凸部分”,凸部分原样保留,非凸部分做一阶展开。这样得到的近似函数可能不再是原问题的严格下界,但仍然是凸函数,可以让子问题有唯一最优解。需要注意的是,如果展开点附近有二阶曲率信息,一阶近似可能偏保守,这时候可以加一个正则项 λ‖x − x(k)‖² 来限制迭代步长,正则项的系数 λ 相当于隐式控制了步长,太大收敛慢,太小又防不住震荡。

4.2 约束凸化:松弛与罚函数

目标函数的近似只是 SCA 的一半,约束的凸化往往更麻烦。非凸等式约束是最难处理的,比如约束 ‖w‖² = 1 这个球面约束,它的集合是非凸的。常见做法是松弛成不等式 ‖w‖² ≤ 1,这一步把约束从等式变不等式后集合变成凸的,但问题来了:松弛后的解可能落在球内部,不满足原等式约束。

应对办法是罚函数法。把等式约束的偏差加入目标函数作为惩罚项,比如加上 μ(‖w‖² − 1)²,然后在外层循环里逐渐增大 μ。这里的 μ 就是罚因子,一开始取 10 或 100,每轮迭代后乘以 1.5,逼迫解逐渐贴近球面。这个技巧在波束成形问题里非常常见,代价是会引入额外的非凸罚项,需要在 SCA 框架内再次把罚项做线性化,相当于两层近似叠加。

另一种常用手段是引入辅助变量。一个看似复杂的非凸约束,往往可以通过新增变量 t 和一个或多个额外约束改写成凸形式。这个技巧在凸优化里叫变量变换,比如把约束 x·y ≥ 1 通过令 u = x, v = y, 约束 u+v 恒定且 uv ≥ 1 来处理,uv ≥ 1 在 u, v 非负时不是凸约束,但通过取对数变成 log u + log v ≥ 0,而 log 函数是凹的,凹函数求和约束大于等于 0 代表的是一个凸集的补集,反而更麻烦。所以变量替换要非常小心,不是所有看起来“换元更简单”的变形都是凸的,每次替换完都要重新验证集合凸性。

4.3 收敛判据与参数选择

SCA 的收敛判断通常检查三样东西:变量变化量、目标函数变化量、梯度或 KKT 残差。变量变化量最常用,即相邻两次迭代满足 ‖x(k+1) − x(k)‖∞ < ε。目标函数变化量适合处理目标值是“只减不增”的单调下降问题,对比前后两轮目标差,小于阈值即停机。KKT 残差是理论最严格的做法,计算量大,工程上见得少。

三个判据对应三种不同场景:变量判据简洁但容易误判,如果目标函数非常平坦,变量还能移动但目标值基本不变了,变量判据会卡住;目标函数判据在非单调场景下容易错过收敛点;混合判据同时检查变量和目标差,虽然代码多两行,但可靠性高得多。参数选择方面,alpha 和 tol 是一对组合,alpha 越小收敛路径越平滑,但最终停留点可能离真最优更远,建议先跑一版 alpha=1 观察曲线是否震荡,如果震荡再降 alpha 到 0.6 左右,同时把 tol 从 1e-6 放宽到 1e-5 看目标值是否变化,如果变化小于 1e-3,精度已经够用。

有一类收敛问题值得单独说:多次随机初始点跑出来的最优值不同,但彼此接近,这属于正常现象,说明问题本身有多个局部驻点。如果最优值相差很大,那就不是收敛判据的问题,而是近似模型构造失配,需要返回 4.1 节检查近似函数。有一点要记住,SCA 并不保证找到全局最优,它给出的是满足 KKT 条件的驻点,实际项目中把驻点配合多次起点挑选出最好的一个,作为工程近似最优解是业界通行的做法。

5. SCA 实现避坑:五条踩坑记录与排查路径

5.1 迭代目标值震荡,曲线呈锯齿状

现象:obj_hist 画出来不是平滑单调下降,而是上下跳动,甚至几步内目标值反复横跳。原因:阻尼系数 alpha 设置过大,或者近似函数在展开点附近与原函数偏差过大,导致子问题最优解越过了原问题可行域的“安全范围”。解决:先把 alpha 降到 0.3 到 0.5 之间,观察两轮曲线是否变平滑;如果仍震荡,检查近似函数是否用了二阶信息而 Hessian 不正定,这种情况要退回一阶展开并增加正则项 λ‖x−x(k)‖²,λ 从 0.1 开始试。

5.2 求解凸子问题时 quadprog 或 CVX 报错,提示问题不可行

现象:在某些迭代轮次,求解器直接返回“Problem is infeasible”,程序中断。原因:这轮里近似约束构造得过于保守,使得约束集合变成了空集。常见于惩罚函数法里 mu 取得过大,或者线性化约束的系数矩阵写错。解决:先打印出错那一轮的约束矩阵,检查是否存在某一行所有系数都为零但右侧常数项大于零;再把罚因子改成渐进式,从 1 开始每轮乘 1.2 而不是直接上大惩罚。

5.3 初始点选得太离谱,收敛到非常差的结果

现象:换一个初始点,最终目标值差了两倍以上,代码逻辑却完全正常。原因:SCA 本质上是局部方法,初始点决定了它落在哪个驻点附近,初始点离可行域太远时,内部子问题的凸近似一开始就偏离了原问题的关键区域。解决:把初始点投影到可行域内再开始迭代,投影可以用 quadprog 先跑一次纯可行解搜索;或者采用“热启动”策略,先用 50 次粗迭代选一组最优中间点,再用这组点作为正式迭代的初始点。

5.4 数值量级差异巨大,迭代到后期精度崩坏

现象:目标函数里同时出现 1e-8 量级和 1e6 量级的项,收敛后目标值明显偏离理论值。原因:MATLAB 默认浮点精度有限,大数吃小数,导致小量级变量在迭代后期几乎不更新。解决:对所有变量和约束做归一化,功率问题里把总功率归一到 1,信道增益归一化到均值 1;如果变量本身跨量级,按列做对角缩放,效果比归一化更细。缩放之后 tol 也要相应调整,原来 1e-6 是绝对量级,归一化后 1e-6 变成了相对精度,判断标准要重新标定。

5.5 收敛了但最终目标值和论文里的结果对不上

现象:曲线收敛得很漂亮,但最终数值和论文或参考实现差了 10% 以上。原因:大概率不是循环写错,而是目标函数本身的表达有差异,比如论文里用的是自然对数 ln,代码里写成了 log10,或者 SINR 分母里加了噪声项而代码漏掉了;还有一种常见情况是约束条件差一个不等号方向,松弛方向和收敛结果完全相反。解决:逐行对照目标函数和约束表达,把论文公式里的每个变量换算关系列成表格,用两三个已知点手动算一遍,对比代码输出。这步没捷径,只能按公式拆开核对。

6. 换到自己问题时怎么验证:多起点检验与目标函数单调性检查

拿到 sca.m 和 xiao_power_beizeng100.m 之后,最常遇到的问题是“我怎么知道我改出来的代码是对的”。我自己的验证套路固定分三步:第一步跑通原示例,确认收敛曲线和目标值和预期一致;第二步把自己的目标函数和约束替换进去,先不追求性能,只验证能跑完整个迭代;第三步做多起点检验,随机生成 20 到 50 个初始点,分别跑 SCA,统计最终目标值的分布。

第三步里有一个可以重复使用的快捷脚本。

% verify_sca.m -- 多起点检验 SCA 稳定性 rng(42); num_starts = 20; final_obj = zeros(num_starts, 1); for i = 1:num_starts x0 = params.xmin + (params.xmax - params.xmin) .* rand(size(params.xmin)); [x_opt, obj_hist] = sca(x0, @approx_func, @solve_func, @obj_func, params); final_obj(i) = obj_func(x_opt); if ~isdecreasing(obj_hist) fprintf('起点 %d: 目标值不单调, 注意震荡风险\n', i); end end fprintf('目标值分布: min=%.6f max=%.6f mean=%.6f\n', ... min(final_obj), max(final_obj), mean(final_obj));

这个脚本里 isdecreasing 是一个自己写的检查函数,判断 obj_hist 序列在最后 20 轮里是否保持单调下降(允许每轮变化小于 1e-4 视为平缓)。做完这一步,如果你的 final_obj 分布集中在很窄的区间里,说明问题结构对初始点不敏感,结果可信度高;如果分布跨度大,说明问题本身多峰严重,或者你的近似函数构造还有问题,需要回去检查近似质量和约束松弛方式。

从那以后我每次把 SCA 代码移植到新问题,都强制先跑一遍这个多起点检验脚本,不通过就不往下继续调参。希望这个验证习惯能帮到你,省掉后面因为结果不可复现而返工的大量时间。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询