MATLAB利用非光滑优化算法解决互补问题

%此函数为主函数，需要输入三个参数，维数n，初始点x0，构造H的辅助向量z，注意x0和z都是列向量；输出x和F。 %函数数据从文件data_g.m，data_y.m，data_z.m获取 %[x F] = ncp(2,[5 5]',[1 1]')在主程序窗口输入这样形式的数据即可。 function [x,F] = ncp(n,x0,z) X = sym(zeros(1,n)); %为符号向量分配空间 max_iter = 5000; % 最大迭代步数 EPS = 1e-6; % 判断算法停止的误差界 %下面是参数 rho = 10; % rho > 0 p= 3; % p > 2 theta = 0.1; % 0 < theta < 1/2 k = 0; xk = x0; % 初始化 for i = 1:1:n X(i) = sym(['x',num2str(i)]); %赋符号向量 end [symF,symgrad_F] = Function(n); %调用函数 Function 得到符号函数symF及其梯度symgrad_F while k < max_iter F = symF; %由于每次迭代xk都会改变，这里需要每次都调用符号函数symF和symgead_F，然后再重新将两函数中的符号X用xk代替 grad_F = symgrad_F; k = k + 1; for j = 1:1:n %用数值代替符号得到函数F和grad_F F = subs(F,X(j),xk(j)); grad_F = subs(grad_F,X(j),xk(j)); end F = eval(F); %将F和grad_F恢复正常精度 grad_F = eval(grad_F); [F_Phi,grad_Phi] = Phi(n,xk,F,grad_F); %调用Phi函数，得到函数F_Phi及其梯度gead_Phi [F_Psi,grad_Psi] = Psi(F_Phi,grad_Phi); %调用Psi函数，得到函数F_Psi及其梯度gead_Psi H = F_H( xk,n,F,grad_F,grad_Phi,z); %调用F_H函数，得到H if norm( grad_Psi ) < EPS %算法停止条件判断 break; else if det( H ) == 0 %H奇异，取负梯度方向 dk = -grad_Psi; else dk = -inv(H)* F_Phi; %H非奇异，像这样取 if grad_Psi' * dk > -rho * (norm(dk))^p %下降条件判断 dk = -grad_Psi; end end ik = get_k( n,xk,dk,symF,F_Psi,grad_Psi,theta); %调用函数get_k，得到最小ik xk = xk + 2^(-ik) * dk; %更新迭代点 end fprintf('This is %dth iteration\n', k); xk end x = xk; %将最后结果赋给x输出 return;