基于逐次凸近似(Successive Convex Approximation)的非凸二次规划问题求解---MATLAB程序

本文主要是介绍基于逐次凸近似(Successive Convex Approximation)的非凸二次规划问题求解---MATLAB程序,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

本文引用了上海财经大学崔雪婷老师最优化理论与方法课程,课程链接如下:

【最优化理论与方法-第十二讲-二次规划】 https://www.bilibili.com/video/BV1vQ4y1P77A/?p=4&share_source=copy_web&vd_source=ec4b99096a4967b6330aae8eaef5e99b

崔老师讲最优化讲的特别好!满分推荐!

逐次凸近似(Successive Convex Approximation, SCA)是一种优化算法,主要应用于求解非凸优化问题。它的基本思想是将一个非凸问题转化为包含多个凸子问题的序列,通过不断的求解凸子问题逼近原问题的最优解。

 图1 非凸函数

现考虑如下非凸二次规划问题,其函数图像如图1所示。

问题1

其中,

原问题的目标函数可以通过特征值分解转化为凸函数减去凸函数的形式,凸函数减去凸函数未必是凸函数

[V,D] = eig(Q);%计算A的特征值对角阵D和特征向量V,使AV=VD成立

其中,矩阵PN都是半正定矩阵,矩阵D的表达式如下所示:

其中,\lambda _{1},\lambda _{2},...,\lambda _{k}\geq 0,\lambda _{k+1},\lambda _{k+2},...< 0

原问题的目标函数可以转化为:

对目标函数的第二项-\left [ x,y \right ]N\left [ x,y \right ]^{T}在点\left ( x^{*},y^{*} \right )处进行凸近似,即在点\left ( x^{*},y^{*} \right )处进行一阶泰勒展开:

至此,原问题可转化为:

 问题2

这样一来,就可以将原来的非凸二次规划问题转化为凸二次规划问题进行求解。

定理:若\left ( x^{*},y^{*} \right )是问题2的最优解,则\left ( x^{*},y^{*} \right )必然是问题1的KKT点(在崔老师的视频中有证明)。

因此,只要找到一个点\left ( x^{*},y^{*} \right )使得\left ( x^{*},y^{*} \right )是问题二的最优解,即可求得原问题的近似最优解。(注意:SCA不能保证得到全局最优解,但解的质量较高)

读到这里,想必各位心中都会有一个疑问:\left ( x^{*},y^{*} \right )点要这么确定呢?SCA算法就是为了找到这样一个点\left ( x^{*},y^{*} \right ),算法步骤如下所示:

1)令k=0,\varepsilon=1\times 10^{-6},取\left ( x_{k},y_{k} \right )\epsilon feasible\: \: region(初值对结果的影响较大,建议取可行域中点);

2)求解近似优化问题(即问题2),得到子问题最优解\left ( x_{k+1},y_{k+1} \right );

3)若\left \| \left ( x_{k+1},y_{k+1} \right )\\ -\, \left ( x_{k},y_{k} \right )\ \right \|\leqslant \varepsilon,输出\left ( x_{k+1},y_{k+1} \right );否则,令k=k+1,转至2)。

MATLAB程序:

clear all
close all
clcQ=[1,0.5;0.5,-1];x=sdpvar(2,1);
xmin=-1;
xmax=1;
Constraints=[];
Constraints=[Constraints,xmin<=x<=xmax];
ops = sdpsettings('solver', 'gurobi', 'verbose', 0);[V,D] = eig(Q);%计算A的特征值对角阵D和特征向量V,使AV=VD成立
lambda_P=D;
lambda_N=-D;
lambda_P(find(D<0))=0;
lambda_N(find(D>0))=0;
P=V*lambda_P*V';
N=V*lambda_N*V';
x0=[0.5;0.5];
x_temp=x0;
while(1)f_k=(x'*P*x-2*x_temp'*N*x+x_temp'*N*x_temp);sol=solvesdp(Constraints,f_k,ops);display([sol.info,' 目标函数值:',num2str(value(x_temp'*Q*x_temp))])x_temp_before=x_temp;x_temp=value(x);if sqrt(sum((x_temp-x_temp_before).^2)/length(x_temp))<1e-10breakend
end
x_result=x_tempX = gridsamp([-1 -1;1 1], 40);
[m,~]=size(X);
YX=zeros(m,1);
for i=1:size(X,1)x=X(i,:);y=x*Q*x';YX(i)=y;
end
X1 = reshape(X(:,1),40,40); X2 = reshape(X(:,2),40,40);
YX = reshape(YX, size(X1));
figure(1), mesh(X1, X2, YX)%绘制预测表面
hold on
scatter3(x_temp(1),x_temp(2),x_temp'*Q*x_temp,200,'r','pentagram','filled')

结果展示:

原问题:

clear all
close all
clcx=sdpvar(2,1);
xmin=-1;
xmax=1;
Constraints=[];
Constraints=[Constraints,xmin<=x<=xmax];
Q=[1,0.5;0.5,-1];
[V,D] = eig(Q);%计算A的特征值对角阵D和特征向量V,使AV=VD成立
lambda_P=D;
lambda_N=-D;
lambda_P(find(D<0))=0;
lambda_N(find(D>0))=0;
P=V*lambda_P*V';
N=V*lambda_N*V';
% P-N
f=x'*Q*x;%x(1)^2-x(2)^2+x(1)*x(2)
% x'*P*x-x'*N*x
ops = sdpsettings('solver', 'gurobi', 'verbose', 0);
sol=solvesdp(Constraints,f,ops);
x=value(x);
display([sol.info,' 目标函数值:',num2str(value(x'*Q*x))])

问题属性:对于非凸二次规划问题,gurobi会将原问题转化为MIP问题进行求解,如图2所示。本文举的例子比较简单,gurobi可以在短时间内求解成功,但对于大规模的非凸二次规划问题,使用gurobi进行求解会面临NP-Hard问题,计算负担较大,用SCA算法可以大大缩短计算时间。

图2

其他子函数:

function  S = gridsamp(range, q)
%GRIDSAMP  n-dimensional grid over given range
%
% Call:    S = gridsamp(range, q)
%
% range :  2*n matrix with lower and upper limits
% q     :  n-vector, q(j) is the number of points
%          in the j'th direction.
%          If q is a scalar, then all q(j) = q
% S     :  m*n array with points, m = prod(q)% hbn@imm.dtu.dk  
% Last update June 25, 2002[mr n] = size(range);    dr = diff(range);
if  mr ~= 2 | any(dr < 0)error('range must be an array with two rows and range(1,:) <= range(2,:)')
end 
sq = size(q);
if  min(sq) > 1 | any(q <= 0)error('q must be a vector with non-negative elements')
end
p = length(q);   
if  p == 1,  q = repmat(q,1,n); 
elseif  p ~= nerror(sprintf('length of q must be either 1 or %d',n))
end % Check for degenerate intervals
i = find(dr == 0);
if  ~isempty(i),  q(i) = 0*q(i); end% Recursive computation
if  n > 1A = gridsamp(range(:,2:end), q(2:end));  % Recursive call[m p] = size(A);   q = q(1);S = [zeros(m*q,1) repmat(A,q,1)];y = linspace(range(1,1),range(2,1), q);k = 1:m;for  i = 1 : qS(k,1) = repmat(y(i),m,1);  k = k + m;end
else    S = linspace(range(1,1),range(2,1), q).';
end

这篇关于基于逐次凸近似(Successive Convex Approximation)的非凸二次规划问题求解---MATLAB程序的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/222788

相关文章

大数据小内存排序问题如何巧妙解决

《大数据小内存排序问题如何巧妙解决》文章介绍了大数据小内存排序的三种方法:数据库排序、分治法和位图法,数据库排序简单但速度慢,对设备要求高;分治法高效但实现复杂;位图法可读性差,但存储空间受限... 目录三种方法:方法概要数据库排序(http://www.chinasem.cn对数据库设备要求较高)分治法(常

Vue项目中Element UI组件未注册的问题原因及解决方法

《Vue项目中ElementUI组件未注册的问题原因及解决方法》在Vue项目中使用ElementUI组件库时,开发者可能会遇到一些常见问题,例如组件未正确注册导致的警告或错误,本文将详细探讨这些问题... 目录引言一、问题背景1.1 错误信息分析1.2 问题原因二、解决方法2.1 全局引入 Element

关于@MapperScan和@ComponentScan的使用问题

《关于@MapperScan和@ComponentScan的使用问题》文章介绍了在使用`@MapperScan`和`@ComponentScan`时可能会遇到的包扫描冲突问题,并提供了解决方法,同时,... 目录@MapperScan和@ComponentScan的使用问题报错如下原因解决办法课外拓展总结@

MybatisGenerator文件生成不出对应文件的问题

《MybatisGenerator文件生成不出对应文件的问题》本文介绍了使用MybatisGenerator生成文件时遇到的问题及解决方法,主要步骤包括检查目标表是否存在、是否能连接到数据库、配置生成... 目录MyBATisGenerator 文件生成不出对应文件先在项目结构里引入“targetProje

C#使用HttpClient进行Post请求出现超时问题的解决及优化

《C#使用HttpClient进行Post请求出现超时问题的解决及优化》最近我的控制台程序发现有时候总是出现请求超时等问题,通常好几分钟最多只有3-4个请求,在使用apipost发现并发10个5分钟也... 目录优化结论单例HttpClient连接池耗尽和并发并发异步最终优化后优化结论我直接上优化结论吧,

Java内存泄漏问题的排查、优化与最佳实践

《Java内存泄漏问题的排查、优化与最佳实践》在Java开发中,内存泄漏是一个常见且令人头疼的问题,内存泄漏指的是程序在运行过程中,已经不再使用的对象没有被及时释放,从而导致内存占用不断增加,最终... 目录引言1. 什么是内存泄漏?常见的内存泄漏情况2. 如何排查 Java 中的内存泄漏?2.1 使用 J

numpy求解线性代数相关问题

《numpy求解线性代数相关问题》本文主要介绍了numpy求解线性代数相关问题,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友们下面随着小编来一起学习学习吧... 在numpy中有numpy.array类型和numpy.mat类型,前者是数组类型,后者是矩阵类型。数组

解决systemctl reload nginx重启Nginx服务报错:Job for nginx.service invalid问题

《解决systemctlreloadnginx重启Nginx服务报错:Jobfornginx.serviceinvalid问题》文章描述了通过`systemctlstatusnginx.se... 目录systemctl reload nginx重启Nginx服务报错:Job for nginx.javas

Redis缓存问题与缓存更新机制详解

《Redis缓存问题与缓存更新机制详解》本文主要介绍了缓存问题及其解决方案,包括缓存穿透、缓存击穿、缓存雪崩等问题的成因以及相应的预防和解决方法,同时,还详细探讨了缓存更新机制,包括不同情况下的缓存更... 目录一、缓存问题1.1 缓存穿透1.1.1 问题来源1.1.2 解决方案1.2 缓存击穿1.2.1

vue解决子组件样式覆盖问题scoped deep

《vue解决子组件样式覆盖问题scopeddeep》文章主要介绍了在Vue项目中处理全局样式和局部样式的方法,包括使用scoped属性和深度选择器(/deep/)来覆盖子组件的样式,作者建议所有组件... 目录前言scoped分析deep分析使用总结所有组件必须加scoped父组件覆盖子组件使用deep前言