利用MATLAB解决非线性规划（NLP）问题

1.利用fmincon函数求解有约束的NLP问题

如果MATLAB中的非线性规划问题可以写成如下形式：

min f(x) Ax<=B A'·x=B' C(x)<=0 C'(x)=0

其中f（x）是标量函数，A，B，A’，B’是相应维数的矩阵和向量，C(x)和C’(x)是非线性向量函数。
x=fmincon(fun,x~0~,a,b,a',b',min,max,nonlcon,options)
它是返回值是向量x，fun是m文件定义的函数 f(x)，a,b,a’,b’是线性约束，如果没有则令a,b,a’,b’=[]；min和max是x的下界和上界，nonlcon是m文件的约束函数，options定义优化参数。

求解下列非线性规划问题：
min f(x)=x₁²+x₂+8
x₁²-x₂>=0
-x₁-x₂²+2=0
x₁,x₂>=0

fun1.m:
function f=fun1(x)
f=x(1)^2+x(2)^2+8;

fun2.m
function [g,h]=fun2(x)
g=-x(1)^2+x(2);
h=-x(1)-x(2)^2+2;

test.m
clc;
options=optimset('LargeScale','off','display','iter');
[x,y]=fmincon(@fun1,[1,1],[],[],[],[],[0,0],[],@fun2,options)

求得x₁=1，x₂=1时，最小值为10.

2.梯度法（最速下降法）

顾名思义，梯度下降法的计算过程就是沿梯度下降的方向求解极小值（也可以沿梯度上升方向求解极大值）。
基本步骤是：对于函数f(x)，给定了一个初始值x0，以及精度ε>0，通过迭代直到 ||∇f(x^k)||²<ε，求步长λ，令x^k+1=x^k+λ∇f(x^k)，再判别此时梯度的范数是否大于ε，求最佳步长的方法：微分法，将带有λ表达式的x^k+1代入到f(x)，令f(x)的导数等于0，求出的λ为最佳步长。
最佳步长表达式一般为(f为一阶偏导数，H为二阶偏导)：

用梯度法求f(x)=4x₁+6x₂-2x₁²-2x₁x₂-2x₂²的极大点，选取初值[1,1]，精度0.01；

方法1，间接法求λ：

fun1.m
function [f,df]=fun1(x)
f=4*x(1)+6*x(2)-2*x(1)^2-2*x(1)*x(2)-2*x(2)^2;
df(1)=4-4*x(1)-2*x(2);
df(2)=6-2*x(1)-4*x(2);

test.m
clc;
x=[1,1];
[f,df]=fun1(x);
while norm(df)^2>0.01syms kx1=x+k*df;[f1,df1]=fun1(x1);s=diff(f1);k=solve(s==0);x=x+k*df[f,df]=fun1(x)
end

求得x= [ 11/32, 21/16 ]，f = 2389/512，

方法2，直接法求λ：

fun1.m
function [f,df,d2f]=fun1(x)
f=4*x(1)+6*x(2)-2*x(1)^2-2*x(1)*x(2)-2*x(2)^2;
df=[4-4*x(1)-2*x(2);6-2*x(1)-4*x(2)];
d2f=[-4 -2;-2 -4];

test.m
clc;
x=[1;1];
[f,df,d2f]=fun1(x);
while norm(df)^2>0.01k=df'*df/(df'*d2f*df)    %k是负的x=x-k*df[f,df,d2f]=fun1(x)
end

方法3，不求λ：

fun1.m
function [f,df]=fun1(x)
f=4*x(1)+6*x(2)-2*x(1)^2-2*x(1)*x(2)-2*x(2)^2;
df=[4-4*x(1)-2*x(2);6-2*x(1)-4*x(2)];

test.m
clc;
x=[1;1];
[f0,df]=fun1(x);
while norm(df)^2>0.01p=df/norm(df);t=1.0;f=fun1(x+t*p);while f<f0t=t/2;f=fun1(x+t*p);endx=x+t*p[f0,df]=fun1(x)
end

3.牛顿法

原理是利用泰勒公式，在x0处展开，且展开到一阶，即f(x) = f(x0)+(x－x0)f’(x0)
求解方程f(x)=0，即f(x0)+(x-x0)f’(x0)=0，求解x = x1=x0－f(x0)/f’(x0)，因为这是利用泰勒公式的一阶展开，f(x) = f(x0)+(x－x0)f’(x0)处并不是完全相等，而是近似相等，这里求得的x1并不能让f（x）=0，只能说f(x1)的值比f(x0)更接近f（x）=0，于是乎，迭代求解的想法就很自然了，可以进而推出x(n+1)=x(n)－f(x(n))/f’(x(n))，通过迭代，这个式子必然在f（x）=0的时候收敛。

在最优化的问题中，线性最优化至少可以使用单纯行法求解，但对于非线性优化问题，牛顿法提供了一种求解的办法。假设任务是优化一个目标函数f，求函数f的极大极小问题，可以转化为求解函数f的导数f’=0的问题，这样求可以把优化问题看成方程求解问题（f’=0）。剩下的问题就和上面提到的牛顿法求解很相似了。

这次为了求解f’=0的根，把f（x）的泰勒展开，展开到2阶形式：

这个式子是成立的，当且仅当 Δx 无线趋近于0。此时上式等价于：

求解：

得出迭代公式：

特别注意：牛顿法对于二次函数极为有效，因为二次函数的二阶偏导矩阵都为常数，一步就可以算出极小值，就不需要选定步长，否则，就按梯度法选择步长的方法。

以2中问题为例

fun1.m
function [f,df,d2f]=fun1(x)
f=4*x(1)+6*x(2)-2*x(1)^2-2*x(1)*x(2)-2*x(2)^2;
df=[4-4*x(1)-2*x(2);6-2*x(1)-4*x(2)];
d2f=[-4 -2;-2 -4];

test.m
clc;
x0=[1;1];
[f,df,d2f]=fun1(x0);
x=x0-d2f\df
f=fun1(x)

考虑非二次函数：
用Newton法求解无约束非线性规划问题：
min f(x)=x₁⁴+25x₂⁴+x₁²x₂²
选取x0=[2;2],ε=10^-6.

fun1.m
function [f,df,d2f]=fun1(x)
f=x(1)^4+25*x(2)^4+x(1)^2*x(2)^2;
df=[4*x(1)^3+2*x(1)*x(2)^2;4*x(1)*x(2)];
d2f=[12*x(1)^2+2*x(2)^2 4*x(1)*x(2);4*x(1)*x(2) 300*x(2)^2+4*x(1)^2];

test.m
clc;
x=[2;2];
[f,df,d2f]=fun1(x);
while norm(df)^2>0.000001k=df'*df/(df'*d2f*df)x=x-k*inv(d2f)*df[f,df,d2f]=fun1(x)
end

注意：如果求极大值，则x=x0+∆；极小值，x=x0-∆