本文主要参考了CSDN上的 另一篇文章，但规范了公式的推导过程和修缮了部分代码

一、代码debug

首先，我们将所有的代码放到MATLAB中，很快在命令行中出现了错误信息

很显然有问题，但是我不知道发生了什么问题。我猜测可能是求解器没有正确安装，因此我正确安装了Gurobi求解器

注意安装Gurobi求解器需要验证license，具体内容可以查询网络上的安装教程。在命令行中输入grbgetkey+licence 就可以完成激活。

但此时还是有问题，我经过查询得知是MEX文件无法指定，是系统路径没有添加gurobi文件的bin，因此我添加到系统路径中

此后文件便可以正确运行了，结果如下。

二、原理

现考虑如下非凸二次规划问题

$$\begin{array}{rlrl}& \min f(x,y)& & =[\begin{array}{c}x,y\end{array}]Q{\left[\begin{array}{c}x,y\end{array}\right]}^T\\ & & & =x^2+xy-y^2\\ & \text{s.t.}& & -1\le x\le 1\\ & & & -1\le y\le 1\end{array}$$

其中，

$$Q=\left[\begin{array}{cc}1& 0.5\\ 0& 5\end{array}\right].$$

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

$$\begin{array}{r}Q=VD{V}^T=V\left({\lambda }_P-{\lambda }_N\right)V^T=\underset{P}{\underbrace{V{\lambda }_P{V}^T}}-\underset{N}{\underbrace{V{\lambda }_N{V}^T}}\end{array}$$

其中，矩阵 $P$ 和 $N$ 都是半正定矩阵，矩阵 $D$ 的表达式如下

$$D=[\begin{array}{cccccc}{\lambda }_1& & & & & \\ & {\lambda }_2& & & & \\ & & \ddots & & & \\ & & & {\lambda }_k& & \\ & & & & {\lambda }_{k+1}& \\ & & & & & \ddots \end{array}]=\underset{{\lambda }_P}{\underbrace{[\begin{array}{cccccc}{\lambda }_1& & & & & \\ & {\lambda }_2& & & & \\ & & \ddots & & & \\ & & & {\lambda }_k& & \\ & & & & 0& \\ & & & & & \ddots \end{array}]}}-\underset{{\lambda }_N}{\underbrace{[\begin{array}{cccccc}0& & & & & \\ & 0& & & & \\ & & \ddots & & & \\ & & & 0& & \\ & & & & -{\lambda }_{k+1}& \\ & & & & & \ddots \end{array}]}}$$

其中 ${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_k\ge 0,{\lambda }_{k+1},{\lambda }_{k+2},\dots <0$。

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

$$\begin{array}{rl}& -\left[x^{\ast },y^{\ast }\right]N{\left[x^{\ast },y^{\ast }\right]}^T-{\left[\nabla \left(\left[x,y\right]N[x,y{]}^T\right)\right|}_{x^{\ast },y^{\ast }}{]}^T{\left(\left[x,y\right]-\left[x^{\ast },y^{\ast }\right]\right)}^T\\ & =-\left[x^{\ast },y^{\ast }\right]N{\left[x^{\ast },y^{\ast }\right]}^T-{\left(2N{\left[x^{\ast },y^{\ast }\right]}^T\right)}^T{\left(\left[x,y\right]-\left[x^{\ast },y^{\ast }\right]\right)}^T\\ & =-2[x^{\ast },y^{\ast }]N[x,y{]}^T+[x^{\ast },y^{\ast }\left]N\right[x^{\ast },y^{\ast }{]}^T\end{array}$$

至此，原问题可转化为：
$$\begin{array}{rlrl}& \min & & f\left(x,y\right)=\left[x,y\right]P{\left[x,y\right]}^T-2\left[x^{\ast },y^{\ast }\right]N{\left[x,y\right]}^T+\left[x^{\ast },y^{\ast }\right]N{\left[x^{\ast },y^{\ast }\right]}^T\\ & s.t.& & -1\le x\le 1\\ & & & -1\le y\le 1\end{array}$$

clear all
close all
clc
 
 
 
Q=[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-10
        break
    end
end
x_result=x_temp
 
X = 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')