非线性规划与matlab解法

Author：怪使
发布时间：August 15, 2019
1853views
No comments
19415 words
Categories：数学建模

非线性规划

1 非线性规划

非线性规划定义：目标函数或者约束条件中包含非线性函数的规划问题。没有通用算法，只能说各个方法都有自己特定的适用范围。

一般形式：

$$ minf(x)\\ {s.t} \ {h_j(x) \le 0,j=1,\cdots,q}\\ \ g_i(x)=0,i=1,\cdots,p\\ $$

其中$x=[{x_1}{\cdots}{x_n}]^T$称为模型的决策变量，$f$称为目标函数，$g_i(i=1,\cdots,p)$和$h_j(j=1,\cdots,q)$称为约束函数。另外的$g_i(x)=0(i=1,\cdots,p)$称为等式约束，${h_j(x) \le 0,j=1,\cdots,q}$称为不等式约束。

线性规划和非线性规划的主要区别：

线性规划如果有最优解，则其最优解只能在可行域的边界上取到；而非线性规划的最优解（若存在），则可能在可行域的任意一点取到。

1.1 非线性规划的Matlab解法：

数学模型写成一般形式为：

$$ minf(x)\\ \begin{cases} Ax \le B\\ Aeq\cdot{x}=Beq\\ C(x) \le 0\\ Ceq(x) =0 \\ \end{cases} $$

其中$f(x)$是标量函数，$A,B,Aeq,Beq$是相应维数的矩阵和向量，$C(x),Ceq(x)$是非线性向量函数。

函数命令是

X=FMINCON(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS)

其返回值是向量$x$,其中FUN是用m文件定义的函数$f(x)$；X0是$x$的初始值，$A,B,Aeq,Beq$定义了线性约束$A*X \le B,Aeq*X=Beq$,如果没有线性约束，可以用[]来替代；LB和UB是变量$x$的下界和上界。NONLCON是用m文件定义的非线性向量函数$C(x),Ceq(x)$；OPTIONS定义了优化参数。

举个例子：求下列非线性规划

$$ minf(x)=x_1^2+x_2^2+x_3^2+8\\ x_1^2-x_2+x_3^2 \ge 0\\ x_1+x_2^2+x_3^3 \le 20\\ -x_1-x_2^2+2=0\\ x_2+2x_3^2=3\\ x_1,x_2,x_3 \ge 0\\ $$

编写m文件fun1.m定义目标函数：

function f=fun1(x);
f = sum(x.^2)+8;

编写m文件fun2.m定义非线性约束条件：

function [g,h]=fun2(x);
g=[-x(1)^2+x(2)-x(3)^2
x(1)+x(2)^2+x(3)^3-20]; %非线性不等式约束
h=[-x(1)-x(2)^2+2
x(2)+2*x(3)^2-3]; %非线性等式约束

编写主程序文件运行：

options=optimset('largescale','off');
[x,y]=fmincon('fun1',rand(3,1),[],[],[],[],zeros(3,1),[], ...
'fun2', options)

得到解为：

$$ x_1=0.5522,x_2=1.2033,x_3= 0.9478,最小值为y=10.6511 $$

1.2 非线性规划问题的迭代求解方法

对于非线性规划模型一般可以采用迭代法进行求最优解:从一个选定的初始点$x^0\in{R^n}$出发，按照某一特定的迭代规则产生一个点列${x_k}$，使得当${x_k}$是有穷点列时，其最后一个点是最优解；当${x_k}$是无穷点列时，它有极限点，其极限点是最优解。

设$x_k\in{R^n}$是某迭代方法的第k轮迭代点，$x_{k+1}\in{R^n}$是第k+1轮迭代点，记

$$ x_{k+1}=x_k+t_kp^k $$

这里$t_k\in{R^1},p^k\in{R^n},\left\|p^k\right\|=1$，并且$p^k$的方向是从点$x_k$向着点$x_{k+1}$的方向。

通常将$p^k$称为第k轮搜索方向，$t_k$为沿$p^k$方向的步长。关键在于如何构造新一轮的搜索方向和确定适当的步长。

设$\overline{x}\in{R^n},p\not= 0$,若存在$\delta \gt 0$,使得：

$$ f(\overline{x}+tp) \lt f(\overline{x}),\forall{t}\in(0,\delta)\\ $$

则称向量$p$为$f$在点$\overline{x}$处的下降方向。

设$\overline{x}\in{R^n},p\not= 0$,若存在$t \gt 0$,使得：

$$ \overline{x}+tp\in{K} $$

称向量$p$为点$\overline{x}$处关于$K$的可行方向。

满足两个条件的，则称为可行下降方向。

基本迭代格式求解非线性规划问题的一般步骤：
选取初始点，构造搜索方向：按照规则构造出可行下降方向，寻求搜索步长：使得目标函数值有某种意义的下降，求出下一个迭代点，直到满足终止条件。

2 无约束问题

2.1 一维搜索方法

当用迭代法求函数的极小点的时候，常常用到一维搜索：即沿某一已知方向求目标函数的极小点。

常见方法：

试探法（斐波那契法）
插值法（抛物线插值法，三次插值法等）
微积分中的求根法（切线法，二分法等）

考虑一维极小化问题

$$ min_{a\le{t}\le{b}}f(t) \tag{1} $$

若$f(t)$是$[a,b]$区间上的下单峰函数，我们通过不断地缩短$[a,b]$的长度，来搜索得到该问题的近似最优解的两个方法。

为了缩短区间$[a,b]$，逐步搜索得到该问题的最优解$t^*$的近似值，可以采用以下途径：在$[a,b]$中任取两个关于$[a,b]$是对称的点$t_1$和$t_2$（不妨设$t_2 \lt t_1$,并把它们叫做搜索点），计算$f(t_1)$和$f(t_2)$并比较它们的大小。对于单峰函数，若$f(t_2) \lt f(t_1)$，则必有$t^*\in[a,t_1]$，因而$[a,t_1]$是缩短了的单峰区间；，若$f(t_1) \lt f(t_2)$，则$[a,t_1]$和$[t_2,b]$都是缩短了的单峰区间。因此通过两个搜索点处目标函数值大小的比较，总可以获得缩短了的单峰区间。对于新的单峰区间重复上述的做法，显然又可以获得更短的单峰区间。如此进行，将单峰区间缩短到充分小的时候，可以取最后的搜索点作为该问题最优解的近似值。

2.1.1 Fibonacci法

具体实现步骤：

选取初始数据，确定单峰区间$[a_0,b_0]$,给出搜索精度$\delta \gt 0$，由$\frac{b-a}{F_n} \le \delta$来确定搜索次数$n$。
$k=1,a=a_0,b=b_0$，计算最初两个搜索点，按照$t_1=a+\frac{F_{n-1}}{F_n}(b-a),t_2=a+\frac{F_{n-2}}{F_n}(b-a)$计算$t_1$和$t_2$。
循环做这件事

$$ while \ k \lt n-1\\ f_1=f(t_1),f_2=f(t_2)\\ if \ f_1 \lt f_2\\ a=t_2;t_2=t_1;t_1=a+\frac{F(n-1-k)}{F(n-k)}(b-a)\\ else\\ b=t_1;t_1=t_2;t_2=a+\frac{F(n-1-k)}{F(n-k)}(a-b) end\\ k=k+1\\ end\\ $$

当进行到$k=n-1$时

$t_1=t_2=\frac{1}{2}(a+b)$

这就无法比较函数值$f(x_1)$和$f(x_2)$的大小而确定最终区间，为此，取到

$$ \begin{cases} t_2=\frac{1}{2}(a+b)\\ t_1=a+(\frac{1}{2}+\varepsilon)(b-a)\\ \end{cases} $$

其中$\varepsilon$为任意小的数。在$t_1$和$t_2$这两点中，以函数值较小者为近似极小点，相应的函数值为近似极小点。并获得最终区间$[a,t_1]$或$[t_2,b]$。

由上述分析可以知道，斐波那契法使用对称搜索的方法，逐步缩短所考察的区间，能够以尽量少的函数求值次数，达到预定的某一缩短率。

2.1.2 0.618法

若$\omega \gt 0$满足比例关系

$$ \frac{\omega}{1}=\frac{1-\omega}{\omega} $$

称之黄金分割数，其值为$\omega=\frac{\sqrt{5}-1}{2}$

黄金分割数和斐波那契分数之间有重要关系：

$$ \omega=\lim_{n \to \infty}\frac{F_{n-1}}{F_n} $$

现在用的是不变的区间缩短率0.618来代替斐波那契法每次不同的缩短率，就可以得到黄金分割法，这个方法可以看作是斐波那契法的近似。

用 0.618 法求解，从第 2 个探索点开始每增加一个探索点作一轮迭代以后，原单峰区间要缩短 0.618 倍。计算 n 个探索点的函数值可以把原区间 $[a_0,b_0]$ 连续缩短 n−1次，因为每次的缩短率均为$\mu$，故最后的区间长度为

$$ (b_0-a_0){\mu}^{n-1} $$

这也就是说，当已知缩短的相对精度为$\delta$时，可以用下面这个式子计算探索点的个数n:

$$ {\mu}^{n-1} \le \delta $$

当然也可以不预先计算探索点的个数n，而是在计算过程中逐次加以判断，看是否满足已给的精度要求。

总结：0.618法是一种等速对称试探法，每次的探索点均取在区间长度的0.618倍和0.382倍处。

2.2 二次插值法

对极小化问题$(1)$，当$f(t)$在$[a,b]$上连续时，可以考虑用多项式插值来进行一维搜索。它的基本思想是：在搜索区间中，不断用低次（通常不超过三次）多项式来近似目标函数，并逐步用插值多项式的极小点来逼近$(1)$的最优解。其实就是缩短搜索区间的拟合函数方法，因为多项式的极小点可比原复杂函数的极小点容易求多了。

2.3 无约束极值问题的解法

无约束极值问题可以表述为

$$ min \quad f(x),x\in{E^n} \tag{2} $$

求解$(2)$的迭代法大体上分为两点：

是用到函数的一阶导数或者二阶导数，称为解析法
是仅仅用到函数值的直接法。

2.3.1 解析法

2.3.1.1 梯度法（最速下降法）

对于基本迭代格式

$$ x_{k+1}=x_k+t_kp^{k} \tag{3} $$

总是考虑从点$x_k$出发沿哪一个方向$p^k$，使目标函数$f$下降得最快。数学分析中提及了：点$x_k$的负梯度方向

$$ p^k=-\nabla{f(x_k)} $$

是从点$x_k$出发使得$f$下降最快的方向，因此，称负梯度方向$-\nabla{f(x_k)}$为$f$在点$x_k$处的最速下降方向。

按照基本迭代格式$(3)$,每一轮从点$x_k$出发沿最速下降方向$-\nabla{f(x_k)}$作一维搜索，来建立求解无约束极值问题的方法，称之为最速下降法。

也就是一个一维搜索方法，然后采用最速梯度下降法来进行。

这个方法的特点是：每轮搜索方向都是目标函数在当前点下降最快的方向。同时，可以用$-\nabla{f(x_k)}=0$或$\left\|-\nabla{f(x_k)}\right\| \le \varepsilon$作为停止条件。

具体步骤如下：

选取初始数据，选取初始点$x_0$,给定终止误差，令$k=0$.
求取梯度向量。计算$\nabla{f(x_k)}$,若$\left\|\nabla{f(x_k)}\right\| \lt \varepsilon$,停止迭代，输出$x_k$。否则，进行3.
构造负梯度方向，取

$$ p^k=-\nabla{f(x_k)} $$

进行一维搜索，求$t_k$，使得

$$ f(x_k+t_kp^k)=min_{t \ge 0}f(x_k+tp^k) $$

令$x_{k+1}=x_k+t_kp^k,k=k+1$,转2.

举个例子：用最速下降法求解无约束非线性规划问题

$$ min \quad f(x)=x_1^2+25x_2^2\\ $$

其中$x=(x_1,x_2)^T$，要求选取初始点$x_0=(2,2)^T$

解：$\nabla{f(x)}=(2x_1,50x_2)^T$

编写M文件detaf.m,定义函数$f(x)$及其梯度列向量如下：

function [f,df]=detaf(x);
f=x(1)^2+25*x(2)^2;  %目标函数
df=[2*x(1)           %梯度函数
50*x(2)];

编写主程序文件如下：

clc
x=[2;2];           %初始点
[f0,g]=detaf(x);  %调用函数。g为梯度函数，f0为目标函数
while norm(g)>0.000001    %误差精细度
 p=-g/norm(g);            %最速下降方向
 t=1.0;f=detaf(x+t*p);       %迭代
 while f>f0             %不满足条件时继续循环
 t=t/2;                    %一维搜索二分
 f=detaf(x+t*p);        %迭代
 end
x=x+t*p;                %迭代
[f0,g]=detaf(x);
end
x,f0

最终得到最小值与最小值的对应点。

2.3.1.2 Newton法

考虑目标函数$f$在点$x_k$处的二次逼近式

$$ f(x)\approx{Q(x)}=f(x_k)+\nabla{f(x_k)^T}(x-x_k)+\frac{1}{2}(x-x_k)^T{\nabla^2}f(x_k)(x-x_k) $$

假定Hesse（海塞）矩阵

$$ \nabla^2{f(x_k)}= \begin{bmatrix} {\frac{{\partial}^2f(x_k)}{\partial{x_1^2}}}&{\cdots}&{\frac{{\partial}^2f(x_k)}{\partial{x_1}{\partial{x_n}}}}\\ {\vdots}&{\ddots}&{\vdots}\\ {\frac{{\partial}^2f(x_k)}{\partial{x_n}{\partial{x_1}}}}&{\cdots}&{\frac{{\partial}^2f(x_k)}{\partial{x_n^2}}}\\ \end{bmatrix} $$

正定

则由于$\nabla^2{f(x_k)}$正定，函数$Q$的驻点$x_{k+1}$是$Q(x)$的极小点。为求这个极小点，我们令：

$$ \nabla{Q(x_{k+1})}=\nabla{f(x_k)}+\nabla^2{f(x_k)}(x_{k+1}-x_k)=0 $$

即可解得到：

$$ x_{k+1}=x_k-[\nabla^2{f(x_k)}]^{-1}\nabla{f(x_k)} $$

对照基本迭代格式（3），可以知道从点$x_k$出发沿着搜索方向

$$ p^k=-[\nabla^2{f(x_k)}]^{-1}\nabla{f(x_k)} $$

并取步长$t_k=1$即可得到$Q(x)$的最小点$x_{k+1}$。通常将方向$p^k$叫做从点$x_k$出发的Newton方向。从一初始点开始，每一轮从当前迭代点出发，沿 Newton 方向并取步长为 1 的求解方法，称之为 Newton 法。

具体步骤如下：

选取初始数据，选取初始点$x_0$,给定终止误差$\varepsilon \gt 0$，令$k=0$.
求取梯度向量。计算$\nabla{f(x_k)}$,若$\left\|\nabla{f(x_k)}\right\| \lt \varepsilon$,停止迭代，输出$x_k$。否则，进行3.
构造Newton方向，计算$[\nabla^2{f(x_k)}]^{-1}$取

$$ p^k=-[\nabla^2{f(x_k)}]^{-1}\nabla{f(x_k)} $$

求下一个迭代点。

令$x_{k+1}=x_k+t_kp^k,k=k+1$,转2.

举个例子：用Newton法求解

$$ min \quad f(x)=x_1^4+25x_2^4+x_1^2x_2^2 $$

选取初始点$x_0=(2,2)^T$

解：

$$ \nabla{f(x)}=[4x_1^3+2x_1x_2^2 \quad 100x_2^3+2x_1^2x_2]^T\\ \nabla^2{f(x)}=\begin{bmatrix} {12x_1^2+2x_2^2}&{4x_1x_2}\\ {4x_1x_2}&{300x_2^2+2x_1^2}\\ \end{bmatrix} $$

编写M文件nwfun.m如下：

function [f,df,d2f]=nwfun(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;100*x(2)^3+2*x(1)^2*x(2)];   %梯度向量
d2f=[2*x(1)^2+2*x(2)^2,4*x(1)*x(2)
 4*x(1)*x(2),300*x(2)^2+2*x(1)^2];   %海塞矩阵向量形式

编写主程序文件并运行：

clc
x=[2;2];               %初始点
[f0,g1,g2]=nwfun(x);   %导入函数
while norm(g1)>0.00001    %精度迭代要求
 p=-inv(g2)*g1;           %迭代，下降方向
 x=x+p;                    %迭代基本格式
 [f0,g1,g2]=nwfun(x);
end
x, f0

如果目标是非二次函数，一般用Newton法通过有限轮迭代并不能保证可以获得最优解。

为了提高计算精度，可以采用可变步长解决这个问题。

clc,clear
x=[2;2];
[f0,g1,g2]=nwfun(x);
while norm(g1)>0.00001
 p=-inv(g2)*g1;p=p/norm(p);   %newton法计算  
 t=1.0;f=nwfun(x+t*p);   
 while f>f0
 t=t/2;f=nwfun(x+t*p);    %更改步长
 end
x=x+t*p;                  %迭代基本格式
[f0,g1,g2]=nwfun(x);
end
x,f0

用该方法的有点就是收敛速度快，缺点则是有时候不太好用。且维数较高时，计算$-[\nabla^2{f(x_k)}]^{-1}$的工作量很大。

2.3.1.3 变尺度法

变尺度法不仅是求解无约束极值问题的有效算法，目前也被推广用来求解有约束的极值问题。因为它既避免了计算二阶导数矩阵（海塞矩阵）及其求逆过程，又比梯度法的收敛速度要快，特别是针对高维问题很有优越性。这里简答介绍一下一种变尺度法——DFP法的基本原理及其计算过程。

我们知道，牛顿法的搜索方向是$-[\nabla^2{f(x_k)}]^{-1}\nabla{f(x_k)}$,为了不计算二阶导数矩阵$[\nabla^2{f(x_k)}]$及其逆矩阵，我们设法构造出一个矩阵来逼近该矩阵的逆矩阵，这类方法也通称为拟牛顿法。

构造这样的近似阵并记为$\overline{H}^{(k)}$。这就要求了：每一步都能以现有的信息来确定下一个搜索方向：每做一次迭代，目标函数值都有所下降，并且这些近似矩阵最后应该收敛于解点出的Hesse矩阵的逆矩阵。

当$f(x)$是二次函数时，其Hesse阵为常数阵$A$，任两点$x_k$和$x_(k+1)$处的梯度之差为

$$ \nabla{f(x_{k+1})}-\nabla{f(x_k)}=A(x_{k+1}-x_k) $$

或

$$ x_{k+1}-x_k=A^{-1}[\nabla{f(x_{k+1})}-\nabla{f(x_k)}] $$

对于非二次函数，仿照二次函数的情形，要求其Hesse矩阵的逆矩阵的第$k+1$次近似矩阵$\overline{H}^{(k+1)}$满足关系式：

$$ x_{k+1}-x_k=\overline{H}^{(k+1)}[\nabla{f(x_{k+1})}-\nabla{f(x_k)}] \tag{4} $$

这就是拟牛顿条件：

如果令

$$ \begin{cases} \Delta{G^{(k)}}=\nabla{f(x_{k+1})}-\nabla{f(x_k)} \\ \Delta{x_k}=x_{k+1}-x_k\\ \end{cases} $$

则式（4）则变为了

$$ \Delta{x_k}=\overline{H}^{(k+1)}\Delta{G^{(k)}} \tag{5} $$

现在假定$\overline{H}^{(k)}$已知，用下面的式子求$\overline{H}^{(k+1)}$(两者皆为对称正定矩阵)

$$ \overline{H}^{(k+1)}=\overline{H}^{(k)}+\Delta\overline{H}^{(k)} $$

其中$\Delta\overline{H}^{(k+1)}$称为第$k$次校正矩阵。显然$\overline{H}^{(k+1)}$应该满足拟牛顿条件（5），也就是要求：

$$ \Delta{x_k}=(\overline{H}^{(k)}+\Delta\overline{H}^{(k)})\Delta{G^{(k)}} $$

或者说：

$$ \Delta\overline{H}^{(k)}\Delta{G^{(k)}} =\Delta{x_k}-\overline{H}^{(k)}\Delta{G^{(k)}} \tag{6} $$

于是可以设想到，$\Delta\overline{H}^{(k)}$的一种比较简单的形式是：

$$ \Delta\overline{H}^{(k)}=\Delta{x_k}(Q^{(k)})^T-\overline{H}^{(k)}\Delta{G^{(k)}}(W^{(k)})^T \tag{7} $$

其中$Q^{k}$与$W^{k}$为两个待定的列向量。

将(7)式子中的$\Delta\overline{H}^{(k)}$代入(6)就可以得到：

$$ \Delta{x_k}(Q^{(k)})^T\Delta{G^{(k)}}-\overline{H}^{(k)}\Delta{G^{(k)}}(W^{(k)})^T \Delta{G^{(k)}} =\Delta{x_k}-\overline{H}^{(k)}\Delta{G^{(k)}} $$

这就说明，应该让：

$$ (Q^{(k)})^T\Delta{G^{(k)}}=(W^{(k)})^T \Delta{G^{(k)}}=1 \tag{8} $$

考虑到$\Delta\overline{H}^{(k)}$应该是对称矩阵，最简单的方式就是取

$$ \begin{cases} Q^{(k)}=\eta_k{\Delta}x^k\\ W^{(k)}=\xi_k\overline{H}^{(k)}\Delta{G^{(k)}}\\ \end{cases} $$

由式子（8）可以得到

$$ \eta_k{({\Delta}x^k)^T}\Delta{G^{(k)}}=\xi_k(\Delta{G^{(k)}})^T\overline{H}^{(k)}\Delta{G^{(k)}}=1 \tag{9} $$

从而求出$\eta_k$与$\xi_k$,然后得到校正矩阵。从而得到尺度矩阵

$$ \overline{H}^{(k+1)}=\overline{H}^{(k)}+\Delta\overline{H}^{(k)} $$

而校正矩阵为：

$$ \Delta\overline{H}^{(k)}=\frac{\Delta{x_k}(\Delta{x_k})^T}{(\Delta{G^{(k)}})^T\Delta{x_k}}-\frac{\overline{H}^{(k)}\Delta{G^{(k)}}(\Delta{G^{(k)}})^T\Delta{H^{(k)}}}{(\Delta{G^{(k)}})^T\overline{H}^{(k)}\Delta{G^{(k)}}} $$

取第一个尺度矩阵$\overline{H}^{(0)}$为单位矩阵，以后的尺度矩阵则按照上面的公式进行逐步形成，可以证明：DFP法的搜索方向为下降方向。

计算步骤：

选取初始数据，选取初始点$x_0$,给定终止误差$\varepsilon \gt 0$，
求取梯度向量。计算$\nabla{f(x_0)}$,若$\left\|\nabla{f(x_0)}\right\| \lt \varepsilon$,则$x_0$为近似极小点，停止迭代，输出$x_k$。否则，进行3.
令

$$ \overline{H}^{(0)}=I(单位矩阵)\\ p^0=-\overline{H}^{(0)}\nabla{f(x_0)}\\ $$

在$p^0$方向进行一维搜索，确定最佳步长$\lambda_0$:

$$ min_\lambda \quad f(x_0+\lambda{p^0})=f(x_0+\lambda_0p^0) $$

如此可以得到下一个近似点

$$ x_1=x_0+\lambda_0p^0 $$

一般地，设已经得到的近似点$x_k$，计算$\nabla{f(x_k)}$,若$\left\|\nabla{f(x_k)}\right\| \lt \varepsilon$,停止迭代，输出$x_k$。否则，计算$\overline{H}^{(k)}$:

$$ \overline{H}^{(k)}=\overline{H}^{(k-1)}+\frac{\Delta{x_{k-1}}(\Delta{x_{k-1}})^T}{(\Delta{G^{(k-1)}})^T\Delta{x_{k-1}}}-\frac{\overline{H}^{(k-1)}\Delta{G^{(k-1)}}(\Delta{G^{(k-1)}})^T\Delta{H^{(k-1)}}}{(\Delta{G^{(k-1)}})^T\overline{H}^{(k-1)}\Delta{G^{(k-1)}}} $$

令$p^k=-\overline{H}^{(k)}\nabla{f(x_k)}$,在$p^k$方向进行一维搜索，确定最佳步长$\lambda_k$,从而得到下一个近似解：

$$ x_{k+1}=x_k+\lambda_kp^k $$

如果$x_{k+1}$满足精度要求，则为所求解，否则回到4，知道求出满足精度的点为止。

2.3.2 直接法

在无约束非线性规划中，遇到的目标函数不可导或者导函数的解析式难以表示的时候，一般使用直接法，这里介绍的是——Powell方法。

该方法主要由：基本搜索、加速搜索、调整搜索方向三部分组成

2.4 Matlab内置函数求无约束极值问题

在内置函数中求解无约束极值问题的函数有fminunc和fminsearch，用法介绍如下：

求函数的极小值

$$ min_x \quad f(x) $$

其中$x$可以是标量或者是向量。

fminunc的基本命令是：

[X,FVAL]=FMINUNC(FUN,X0,OPTIONS,P1,P2,...)

返回函数值X是极小点，FVAL是函数的极小值。FUN是一个m文件，当FUN只有一个返回值的时候，返回值是目标函数$f(x)$，当有两个返回值时，第二个返回值是$f(x)$的梯度向量；当有三个返回值的时候，第三个返回值是Hessian矩阵。X0是向量$x$的初始值。OPTIONS是优化参数，而P1,P2是可以传递给FUN的一些参数（含参数的函数之类的）

举个例子：求函数$f(x)=100(x_2-x_1^2)^2+(1-x_1)^2$的最小值、

编写m文件fun2.m，如下

function [f,g]=fun2(x);
f=100*(x(2)-x(1)^2)^2+(1-x(1))^2;  %目标函数
g=[-400*x(1)*(x(2)-x(1)^2)-2*(1-x(1));200*(x(2)-x(1)^2)];   %梯度函数

编写主函数文件并运行：

options = optimset('GradObj' , 'on' );  %优化参数
[x,y]=fminunc('fun2' ,rand(1,2),options)

即可求得函数的极小值，当然你也可以利用二阶导数：

function [f,df,d2f]=fun3(x);
f=100*(x(2)-x(1)^2)^2+(1-x(1))^2;
df=[-400*x(1)*(x(2)-x(1)^2)-2*(1-x(1));200*(x(2)-x(1)^2)];
d2f=[-400*x(2)+1200*x(1)^2+2,-400*x(1)
 -400*x(1),200];

然后主程序的优化参数记得打开海塞矩阵：

options = optimset('GradObj' , 'on' , 'Hessian' , 'on' );
[x,y]=fminunc('fun3' ,rand(1,2),options)

求多元函数的极值也可以使用fminsearch命令，其使用格式为：

[X,FVAL,EXITFLAG,IUTPUT]=FMINSEARCH(FUN,X0,OPTIONS,P1,P2,...)

例如：求函数$f(x)=sin(x)+3$取到最小值的$x$值

编写M文件fun4.m如下：

function f=fun4(x);
f=sin(x)+3;

主程序如下：

x0=2;
[x,y]=fminsearch(@fun4,x0)

即可求得在初值2附近的极小点和极小值。

3 带约束的极值问题

带有约束条件的极值问题称为约束极值问题，也叫做规划问题。

思路：将约束问题转化为无约束问题，将非线性规划转化为线性规划，主要就是把复杂变为简单。

库恩——塔克条件是非线性规划的重点，是确定某点是否为最优点的必要条件，一般来说不是充分条件（除了凸规划）

3.1 二次规划

若某非线性规划的目标函数为自变量$x$的二次函数，约束条件又全部是线性的，就称为二次规划

数学模型：

$$ min \quad \frac{1}{2}x^THx+f^Tx,\\ s.t.\quad \begin{cases} Ax \le b\\ Aeq{\cdot}x=beq\\ \end{cases} $$

这里$H$是是对称矩阵，$f,b$是列向量，$A$是相应维数的矩阵。

求解二次规划的命令是：

[X,FVAL]=QUADPROG(H,f,A,b,Aeq,beq,LB,UB,X0,OPTIONS)

返回值X是决策向量$x$的值。

举个例子：求解二次规划

$$ \begin{cases} min \quad f(x)=2x_1^2-4x_1x_2+4x_2^2-6x_1-3x_2\\ x_1+x_2 \le 3\\ 4x_1+x_2 \le 9\\ x_1,x_2 \ge 0\\ \end{cases} $$

编程如下：

h=[4,-4;-4,8];
f=[-6;-3];
a=[1,1;4,1];
b=[3;9];
[x,value]=quadprog(h,f,a,b,[],[],zeros(2,1))  %下界大于等于0，这里给了不等式约束，等式约束没有给了[]

3.2 罚函数法

利用罚函数法，可以将非线性规划的问题转化为无约束极值问题求解，这种方法被称为序列无约束最小化技术，简称SUMT。

主要思想：利用问题中的约束函数做出适当的罚函数，由此构造出带有参数的增广目标函数，从而把问题转化为无约束非线性规划问题。两种：内与外罚函数法。

外罚函数法：

$$ min \quad f(x)\\ s.t.\begin{cases} g_i(x) \le 0,\quad i=1,\cdots,r,\\ h_j(x) \ge 0,\quad j=1,\cdots,s,\\ k_m(x)=0,\quad m=1,\cdots,t\\ \end{cases} $$

取一个充分大的数$M \gt 0$，构造函数

$$ P(x,M)=f(x)+M\sum^r_{i=1}max(g_i(x),0)-M\sum^s_{i=1}min(h_i(x),0)+M\sum^t_{i=1}|k_i(x)| $$

以增广目标函数$P(x,M)$为目标函数的无约束极值问题

$$ min \quad P(x,M) $$

的最优解$x$也是原问题的最优解。

举个例子：求下列非线性规划

$$ \begin{cases} min \ f(x)=x_1^2+x_2^2+8\\ x_1^2-x_2 \ge 0\\ -x_1-x_2^2+2=0 \\ x_1,x_2 \ge 0\\ \end{cases} $$

解：编写M文件test.m

function g=test(x);
M=50000;   %M取到足够大
f=x(1)^2+x(2)^2+8;
g=f-M*min(x(1),0)-M*min(x(2),0)-M*min(x(1)^2-x(2),0)+...
    M*abs(-x(1)-x(2)^2+2);

在窗口运行

 [x,y]=fminunc('test',rand(2,1))

即可获得问题的解。

3.3 Matlab内置函数求解约束极值问题

这些内置函数有很多：fminbnd、fmincon、quadprog、fseminf、fminimax，之前已经介绍了函数fmincon和quadprog

3.3.1 fminbnd函数

求单变量非线性函数在区间上的极小值

$$ min_xf(x),x\in [x_1,x_2] $$

命令为：

[X,FVAL]=FMINBND(FUN,X1,X2,OPTIONS)

返回值是极小点$x$与函数的极小值，这里的FUN是m文件定义的函数或者是内置的单变量数学函数。

举个例子：求函数$f(x)=(x-3)^2-1,x\in[0,5]$的最小值

编写m文件fun5.m

function f=fun5(x);
 f=(x-3)^2-1;

在窗口输入

[x,y]=fminbnd('fun5',0,5)

即可求得极小点和极小值

3.3.2 fseminf函数

求

$$ min_x\left\{ F(x)|C(x) \le 0,Ceq(x)=0,PHI(x,w) \le 0 \right\}\\ s.t. \begin{cases} A*x \le B\\ Aeq*x=Beq\\ \end{cases} $$

其中$C(x),Ceq(x),PHI(x,w)$都是向量函数，$w$是附加的向量变量，$w$的每个分量都限定在某个区间内。

命令格式为：

X=FSEMINF(FUN,X0,NTHETA,SEMINFCON,A,B,Aeq，Beq)

其中FUN用于定义目标函数$f(x)$,X0为$x$的初始值；NTHETA是半无穷约束PHI$(x,w)$的个数；函数SEMINFCON用于定义非线性不等式约束$C(x)$,非线性等式约束$Ceq(x)$和半无穷约束$PHI(x,w)$的每一个分量函数，函数SEMINFCON有两个输入参量$X$和$S$,$S$是推荐的取样步长，可以取用。

举个例子：求函数$f(x)=(x_1-0.5)^2+(x_2-0.5)^2+(x_3-0.5)^2$取最小值时的$x$值，约束为：

$$ K_1(x,w_1)=sin(w_1,x_1)cos(w_1,x_2)-\frac{1}{1000}(w_1-50)^2-sin(w_1x_3)-x_3 \le 1\\ K_2(x,w_2)=sin(w_2,x_2)cos(w_2,x_1)-\frac{1}{1000}(w_2-50)^2-sin(w_2x_3)-x_3 \le 1\\ 1 \le w_1 \le 100, \ 1 \le w_2 \le 100\\ $$

编写M文件fun6.m定义目标函数如下：

function f=fun6(x,s);
f=sum((x-0.5).^2);

编写M文件fun7.m定义约束条件如下：

function [c,ceq,k1,k2,s]=fun7(x,s);
c=[];ceq=[];
if isnan(s(1,1))
 s=[0.2,0;0.2 0];
end
%取样值
w1=1:s(1,1):100;
w2=1:s(2,1):100;
%半无穷约束
k1=sin(w1*x(1)).*cos(w1*x(2))-1/1000*(w1-50).^2-sin(w1*x(3))-x(3)-1;
k2=sin(w2*x(2)).*cos(w2*x(1))-1/1000*(w2-50).^2-sin(w2*x(3))-x(3)-1;
%画出半无穷约束的图形
plot(w1,k1,'-',w2,k2,'+');

在命令行窗口键入：

[x,y]=fseminf(@fun6,rand(3,1),2,@fun7)  %给与初始点和约束个数

3.3.3 fminimax函数

求解

$$ min_x\left\{max_{F_i}F(x)\right\}\\ s.t.\begin{cases} A*x \le b\\ Aeq*x=Beq \\ C(x) \le 0\\ Ceq(x) = 0 \\ LB \le x \le UB\\ \end{cases} $$

其中$F(x)=\left\{F_1(x),\cdots,F_m(x)\right\}$

上述问题的命令为：

X=FMINIMAX(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON)

举个例子：求函数族$\left\{f_1(x),f_2(x),f_3(x),f_4(x),f_5(x)\right\}$取极大极小值时的$x$值。其中：

$$ \begin{cases} f_1(x)=2x_1^2+x_2^2-48x_1-40x_2+304\\ f_2(x)=-x_1^2-3x_2^2\\ f_3(x)=x_1+3x_2-18\\ f_4(x)=-x_1-x_2\\ f_5(x)=x_1+x_2-8\\ \end{cases} $$

编写M文件fun8.m定义向量函数如下：

function f=fun8(x);
f=[2*x(1)^2+x(2)^2-48*x(1)-40*x(2)+304
 -x(1)^2-3*x(2)^2
 x(1)+3*x(2)-18
 -x(1)-x(2)
 x(1)+x(2)-8];

调用函数

[x,y]=fminimax(@fun8,rand(2,1))

3.3.4 利用梯度求解约束优化问题

举个例子：已知函数$f(x)=e^{x_1}(4x_1^2+2x_2^2+4x_1x_2+2x_2+1)$,且满足非线性约束：

$$ \begin{cases} x_1x_2-x_1-x_2 \le -1.5\\ x_1x_2 \gt -10\\ \end{cases} $$

求$min_xf(x)$

题目中的目标函数梯度为：

$$ \begin{bmatrix} e^{x_1}(4x_1^2+2x_2^2+4x_1x_2+8x_1+6x_2+1)\\ e^{x_1}(4x_1+4x_2+2)\\ \end{bmatrix} $$

编写M文件fun9.m定义目标函数及梯度函数：

function [f,df]=fun9(x)
f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);
df=[exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+8*x(1)+6*x(2)+1);exp(x(1))*(4*x(2)+4*x(1)+2)];

编写M文件fun10.m定义约束条件及约束条件的梯度函数

function [c,ceq,dc,dceq]=fun10(x)
c=[x(1)*x(2)-x(1)-x(2)+1.5;-x(1)*x(2)-10];
dc=[x(2)-1,-x(2);x(1)-1,-x(1)];
ceq=[];dceq=[];

调用函数fmincon，编写主函数并运行：

%采用标准算法
options=optimset('largescale','off');
%采用梯度
options=optimset(options,'GradObj','on','GradConstr','on');
[x,y]=fmincon(@fun9,rand(2,1),[],[],[],[],[],[],@fun10,options)

3.4 优化工具箱的用户图形界面解法

Matlab 优化工具箱中的 optimtool 命令提供了优化问题的用户图形界面解法。optimtool 可应用到所有优化问题的求解，计算结果可以输出到 Matlab 工作空间中。

Last modification：August 15th, 2019 at 05:48 pm

如果觉得我的文章对你有用，请随意赞赏

非线性规划与matlab解法

怪使 • 2019 年 08 月 15 日

非线性规划

1 非线性规划

非线性规划定义：目标函数或者约束条件中包含非线性函数的规划问题。没有通用算法，只能说各个方法都有自己特定的适用范围。

一般形式：

$$ minf(x)\\ {s.t} \ {h_j(x) \le 0,j=1,\cdots,q}\\ \ g_i(x)=0,i=1,\cdots,p\\ $$

线性规划和非线性规划的主要区别：

线性规划如果有最优解，则其最优解只能在可行域的边界上取到；而非线性规划的最优解（若存在），则可能在可行域的任意一点取到。

1.1 非线性规划的Matlab解法：

数学模型写成一般形式为：

$$ minf(x)\\ \begin{cases} Ax \le B\\ Aeq\cdot{x}=Beq\\ C(x) \le 0\\ Ceq(x) =0 \\ \end{cases} $$

其中$f(x)$是标量函数，$A,B,Aeq,Beq$是相应维数的矩阵和向量，$C(x),Ceq(x)$是非线性向量函数。

函数命令是

X=FMINCON(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS)

举个例子：求下列非线性规划

$$ minf(x)=x_1^2+x_2^2+x_3^2+8\\ x_1^2-x_2+x_3^2 \ge 0\\ x_1+x_2^2+x_3^3 \le 20\\ -x_1-x_2^2+2=0\\ x_2+2x_3^2=3\\ x_1,x_2,x_3 \ge 0\\ $$

编写m文件fun1.m定义目标函数：

function f=fun1(x);
f = sum(x.^2)+8;

编写m文件fun2.m定义非线性约束条件：

function [g,h]=fun2(x);
g=[-x(1)^2+x(2)-x(3)^2
x(1)+x(2)^2+x(3)^3-20]; %非线性不等式约束
h=[-x(1)-x(2)^2+2
x(2)+2*x(3)^2-3]; %非线性等式约束

编写主程序文件运行：

options=optimset('largescale','off');
[x,y]=fmincon('fun1',rand(3,1),[],[],[],[],zeros(3,1),[], ...
'fun2', options)

得到解为：

$$ x_1=0.5522,x_2=1.2033,x_3= 0.9478,最小值为y=10.6511 $$

1.2 非线性规划问题的迭代求解方法

设$x_k\in{R^n}$是某迭代方法的第k轮迭代点，$x_{k+1}\in{R^n}$是第k+1轮迭代点，记

$$ x_{k+1}=x_k+t_kp^k $$

这里$t_k\in{R^1},p^k\in{R^n},\left\|p^k\right\|=1$，并且$p^k$的方向是从点$x_k$向着点$x_{k+1}$的方向。

通常将$p^k$称为第k轮搜索方向，$t_k$为沿$p^k$方向的步长。关键在于如何构造新一轮的搜索方向和确定适当的步长。

设$\overline{x}\in{R^n},p\not= 0$,若存在$\delta \gt 0$,使得：

$$ f(\overline{x}+tp) \lt f(\overline{x}),\forall{t}\in(0,\delta)\\ $$

则称向量$p$为$f$在点$\overline{x}$处的下降方向。

设$\overline{x}\in{R^n},p\not= 0$,若存在$t \gt 0$,使得：

$$ \overline{x}+tp\in{K} $$

称向量$p$为点$\overline{x}$处关于$K$的可行方向。

满足两个条件的，则称为可行下降方向。

2 无约束问题

2.1 一维搜索方法

当用迭代法求函数的极小点的时候，常常用到一维搜索：即沿某一已知方向求目标函数的极小点。

常见方法：

试探法（斐波那契法）
插值法（抛物线插值法，三次插值法等）
微积分中的求根法（切线法，二分法等）

考虑一维极小化问题

$$ min_{a\le{t}\le{b}}f(t) \tag{1} $$

若$f(t)$是$[a,b]$区间上的下单峰函数，我们通过不断地缩短$[a,b]$的长度，来搜索得到该问题的近似最优解的两个方法。

2.1.1 Fibonacci法

具体实现步骤：

选取初始数据，确定单峰区间$[a_0,b_0]$,给出搜索精度$\delta \gt 0$，由$\frac{b-a}{F_n} \le \delta$来确定搜索次数$n$。
$k=1,a=a_0,b=b_0$，计算最初两个搜索点，按照$t_1=a+\frac{F_{n-1}}{F_n}(b-a),t_2=a+\frac{F_{n-2}}{F_n}(b-a)$计算$t_1$和$t_2$。
循环做这件事

$$ while \ k \lt n-1\\ f_1=f(t_1),f_2=f(t_2)\\ if \ f_1 \lt f_2\\ a=t_2;t_2=t_1;t_1=a+\frac{F(n-1-k)}{F(n-k)}(b-a)\\ else\\ b=t_1;t_1=t_2;t_2=a+\frac{F(n-1-k)}{F(n-k)}(a-b) end\\ k=k+1\\ end\\ $$

当进行到$k=n-1$时

$t_1=t_2=\frac{1}{2}(a+b)$

这就无法比较函数值$f(x_1)$和$f(x_2)$的大小而确定最终区间，为此，取到

$$ \begin{cases} t_2=\frac{1}{2}(a+b)\\ t_1=a+(\frac{1}{2}+\varepsilon)(b-a)\\ \end{cases} $$

其中$\varepsilon$为任意小的数。在$t_1$和$t_2$这两点中，以函数值较小者为近似极小点，相应的函数值为近似极小点。并获得最终区间$[a,t_1]$或$[t_2,b]$。

由上述分析可以知道，斐波那契法使用对称搜索的方法，逐步缩短所考察的区间，能够以尽量少的函数求值次数，达到预定的某一缩短率。

2.1.2 0.618法

若$\omega \gt 0$满足比例关系

$$ \frac{\omega}{1}=\frac{1-\omega}{\omega} $$

称之黄金分割数，其值为$\omega=\frac{\sqrt{5}-1}{2}$

黄金分割数和斐波那契分数之间有重要关系：

$$ \omega=\lim_{n \to \infty}\frac{F_{n-1}}{F_n} $$

现在用的是不变的区间缩短率0.618来代替斐波那契法每次不同的缩短率，就可以得到黄金分割法，这个方法可以看作是斐波那契法的近似。

$$ (b_0-a_0){\mu}^{n-1} $$

这也就是说，当已知缩短的相对精度为$\delta$时，可以用下面这个式子计算探索点的个数n:

$$ {\mu}^{n-1} \le \delta $$

当然也可以不预先计算探索点的个数n，而是在计算过程中逐次加以判断，看是否满足已给的精度要求。

总结：0.618法是一种等速对称试探法，每次的探索点均取在区间长度的0.618倍和0.382倍处。

2.2 二次插值法

2.3 无约束极值问题的解法

无约束极值问题可以表述为

$$ min \quad f(x),x\in{E^n} \tag{2} $$

求解$(2)$的迭代法大体上分为两点：

是用到函数的一阶导数或者二阶导数，称为解析法
是仅仅用到函数值的直接法。

2.3.1 解析法

2.3.1.1 梯度法（最速下降法）

对于基本迭代格式

$$ x_{k+1}=x_k+t_kp^{k} \tag{3} $$

总是考虑从点$x_k$出发沿哪一个方向$p^k$，使目标函数$f$下降得最快。数学分析中提及了：点$x_k$的负梯度方向

$$ p^k=-\nabla{f(x_k)} $$

是从点$x_k$出发使得$f$下降最快的方向，因此，称负梯度方向$-\nabla{f(x_k)}$为$f$在点$x_k$处的最速下降方向。

按照基本迭代格式$(3)$,每一轮从点$x_k$出发沿最速下降方向$-\nabla{f(x_k)}$作一维搜索，来建立求解无约束极值问题的方法，称之为最速下降法。

也就是一个一维搜索方法，然后采用最速梯度下降法来进行。

具体步骤如下：

选取初始数据，选取初始点$x_0$,给定终止误差，令$k=0$.
求取梯度向量。计算$\nabla{f(x_k)}$,若$\left\|\nabla{f(x_k)}\right\| \lt \varepsilon$,停止迭代，输出$x_k$。否则，进行3.
构造负梯度方向，取

$$ p^k=-\nabla{f(x_k)} $$

进行一维搜索，求$t_k$，使得

$$ f(x_k+t_kp^k)=min_{t \ge 0}f(x_k+tp^k) $$

令$x_{k+1}=x_k+t_kp^k,k=k+1$,转2.

举个例子：用最速下降法求解无约束非线性规划问题

$$ min \quad f(x)=x_1^2+25x_2^2\\ $$

其中$x=(x_1,x_2)^T$，要求选取初始点$x_0=(2,2)^T$

解：$\nabla{f(x)}=(2x_1,50x_2)^T$

编写M文件detaf.m,定义函数$f(x)$及其梯度列向量如下：

function [f,df]=detaf(x);
f=x(1)^2+25*x(2)^2;  %目标函数
df=[2*x(1)           %梯度函数
50*x(2)];

编写主程序文件如下：

clc
x=[2;2];           %初始点
[f0,g]=detaf(x);  %调用函数。g为梯度函数，f0为目标函数
while norm(g)>0.000001    %误差精细度
 p=-g/norm(g);            %最速下降方向
 t=1.0;f=detaf(x+t*p);       %迭代
 while f>f0             %不满足条件时继续循环
 t=t/2;                    %一维搜索二分
 f=detaf(x+t*p);        %迭代
 end
x=x+t*p;                %迭代
[f0,g]=detaf(x);
end
x,f0

最终得到最小值与最小值的对应点。

2.3.1.2 Newton法

考虑目标函数$f$在点$x_k$处的二次逼近式

$$ f(x)\approx{Q(x)}=f(x_k)+\nabla{f(x_k)^T}(x-x_k)+\frac{1}{2}(x-x_k)^T{\nabla^2}f(x_k)(x-x_k) $$

假定Hesse（海塞）矩阵

正定

则由于$\nabla^2{f(x_k)}$正定，函数$Q$的驻点$x_{k+1}$是$Q(x)$的极小点。为求这个极小点，我们令：

$$ \nabla{Q(x_{k+1})}=\nabla{f(x_k)}+\nabla^2{f(x_k)}(x_{k+1}-x_k)=0 $$

即可解得到：

$$ x_{k+1}=x_k-[\nabla^2{f(x_k)}]^{-1}\nabla{f(x_k)} $$

对照基本迭代格式（3），可以知道从点$x_k$出发沿着搜索方向

$$ p^k=-[\nabla^2{f(x_k)}]^{-1}\nabla{f(x_k)} $$

具体步骤如下：

选取初始数据，选取初始点$x_0$,给定终止误差$\varepsilon \gt 0$，令$k=0$.
求取梯度向量。计算$\nabla{f(x_k)}$,若$\left\|\nabla{f(x_k)}\right\| \lt \varepsilon$,停止迭代，输出$x_k$。否则，进行3.
构造Newton方向，计算$[\nabla^2{f(x_k)}]^{-1}$取

$$ p^k=-[\nabla^2{f(x_k)}]^{-1}\nabla{f(x_k)} $$

求下一个迭代点。

令$x_{k+1}=x_k+t_kp^k,k=k+1$,转2.

举个例子：用Newton法求解

$$ min \quad f(x)=x_1^4+25x_2^4+x_1^2x_2^2 $$

选取初始点$x_0=(2,2)^T$

解：

$$ \nabla{f(x)}=[4x_1^3+2x_1x_2^2 \quad 100x_2^3+2x_1^2x_2]^T\\ \nabla^2{f(x)}=\begin{bmatrix} {12x_1^2+2x_2^2}&{4x_1x_2}\\ {4x_1x_2}&{300x_2^2+2x_1^2}\\ \end{bmatrix} $$

编写M文件nwfun.m如下：

function [f,df,d2f]=nwfun(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;100*x(2)^3+2*x(1)^2*x(2)];   %梯度向量
d2f=[2*x(1)^2+2*x(2)^2,4*x(1)*x(2)
 4*x(1)*x(2),300*x(2)^2+2*x(1)^2];   %海塞矩阵向量形式

编写主程序文件并运行：

clc
x=[2;2];               %初始点
[f0,g1,g2]=nwfun(x);   %导入函数
while norm(g1)>0.00001    %精度迭代要求
 p=-inv(g2)*g1;           %迭代，下降方向
 x=x+p;                    %迭代基本格式
 [f0,g1,g2]=nwfun(x);
end
x, f0

如果目标是非二次函数，一般用Newton法通过有限轮迭代并不能保证可以获得最优解。

为了提高计算精度，可以采用可变步长解决这个问题。

clc,clear
x=[2;2];
[f0,g1,g2]=nwfun(x);
while norm(g1)>0.00001
 p=-inv(g2)*g1;p=p/norm(p);   %newton法计算  
 t=1.0;f=nwfun(x+t*p);   
 while f>f0
 t=t/2;f=nwfun(x+t*p);    %更改步长
 end
x=x+t*p;                  %迭代基本格式
[f0,g1,g2]=nwfun(x);
end
x,f0

用该方法的有点就是收敛速度快，缺点则是有时候不太好用。且维数较高时，计算$-[\nabla^2{f(x_k)}]^{-1}$的工作量很大。

2.3.1.3 变尺度法

当$f(x)$是二次函数时，其Hesse阵为常数阵$A$，任两点$x_k$和$x_(k+1)$处的梯度之差为

$$ \nabla{f(x_{k+1})}-\nabla{f(x_k)}=A(x_{k+1}-x_k) $$

或

$$ x_{k+1}-x_k=A^{-1}[\nabla{f(x_{k+1})}-\nabla{f(x_k)}] $$

对于非二次函数，仿照二次函数的情形，要求其Hesse矩阵的逆矩阵的第$k+1$次近似矩阵$\overline{H}^{(k+1)}$满足关系式：

$$ x_{k+1}-x_k=\overline{H}^{(k+1)}[\nabla{f(x_{k+1})}-\nabla{f(x_k)}] \tag{4} $$

这就是拟牛顿条件：

如果令

$$ \begin{cases} \Delta{G^{(k)}}=\nabla{f(x_{k+1})}-\nabla{f(x_k)} \\ \Delta{x_k}=x_{k+1}-x_k\\ \end{cases} $$

则式（4）则变为了

$$ \Delta{x_k}=\overline{H}^{(k+1)}\Delta{G^{(k)}} \tag{5} $$

现在假定$\overline{H}^{(k)}$已知，用下面的式子求$\overline{H}^{(k+1)}$(两者皆为对称正定矩阵)

$$ \overline{H}^{(k+1)}=\overline{H}^{(k)}+\Delta\overline{H}^{(k)} $$

其中$\Delta\overline{H}^{(k+1)}$称为第$k$次校正矩阵。显然$\overline{H}^{(k+1)}$应该满足拟牛顿条件（5），也就是要求：

$$ \Delta{x_k}=(\overline{H}^{(k)}+\Delta\overline{H}^{(k)})\Delta{G^{(k)}} $$

或者说：

$$ \Delta\overline{H}^{(k)}\Delta{G^{(k)}} =\Delta{x_k}-\overline{H}^{(k)}\Delta{G^{(k)}} \tag{6} $$

于是可以设想到，$\Delta\overline{H}^{(k)}$的一种比较简单的形式是：

$$ \Delta\overline{H}^{(k)}=\Delta{x_k}(Q^{(k)})^T-\overline{H}^{(k)}\Delta{G^{(k)}}(W^{(k)})^T \tag{7} $$

其中$Q^{k}$与$W^{k}$为两个待定的列向量。

将(7)式子中的$\Delta\overline{H}^{(k)}$代入(6)就可以得到：

$$ \Delta{x_k}(Q^{(k)})^T\Delta{G^{(k)}}-\overline{H}^{(k)}\Delta{G^{(k)}}(W^{(k)})^T \Delta{G^{(k)}} =\Delta{x_k}-\overline{H}^{(k)}\Delta{G^{(k)}} $$

这就说明，应该让：

$$ (Q^{(k)})^T\Delta{G^{(k)}}=(W^{(k)})^T \Delta{G^{(k)}}=1 \tag{8} $$

考虑到$\Delta\overline{H}^{(k)}$应该是对称矩阵，最简单的方式就是取

$$ \begin{cases} Q^{(k)}=\eta_k{\Delta}x^k\\ W^{(k)}=\xi_k\overline{H}^{(k)}\Delta{G^{(k)}}\\ \end{cases} $$

由式子（8）可以得到

$$ \eta_k{({\Delta}x^k)^T}\Delta{G^{(k)}}=\xi_k(\Delta{G^{(k)}})^T\overline{H}^{(k)}\Delta{G^{(k)}}=1 \tag{9} $$

从而求出$\eta_k$与$\xi_k$,然后得到校正矩阵。从而得到尺度矩阵

$$ \overline{H}^{(k+1)}=\overline{H}^{(k)}+\Delta\overline{H}^{(k)} $$

而校正矩阵为：

取第一个尺度矩阵$\overline{H}^{(0)}$为单位矩阵，以后的尺度矩阵则按照上面的公式进行逐步形成，可以证明：DFP法的搜索方向为下降方向。

计算步骤：

选取初始数据，选取初始点$x_0$,给定终止误差$\varepsilon \gt 0$，
求取梯度向量。计算$\nabla{f(x_0)}$,若$\left\|\nabla{f(x_0)}\right\| \lt \varepsilon$,则$x_0$为近似极小点，停止迭代，输出$x_k$。否则，进行3.
令

$$ \overline{H}^{(0)}=I(单位矩阵)\\ p^0=-\overline{H}^{(0)}\nabla{f(x_0)}\\ $$

在$p^0$方向进行一维搜索，确定最佳步长$\lambda_0$:

$$ min_\lambda \quad f(x_0+\lambda{p^0})=f(x_0+\lambda_0p^0) $$

如此可以得到下一个近似点

$$ x_1=x_0+\lambda_0p^0 $$

一般地，设已经得到的近似点$x_k$，计算$\nabla{f(x_k)}$,若$\left\|\nabla{f(x_k)}\right\| \lt \varepsilon$,停止迭代，输出$x_k$。否则，计算$\overline{H}^{(k)}$:

令$p^k=-\overline{H}^{(k)}\nabla{f(x_k)}$,在$p^k$方向进行一维搜索，确定最佳步长$\lambda_k$,从而得到下一个近似解：

$$ x_{k+1}=x_k+\lambda_kp^k $$

如果$x_{k+1}$满足精度要求，则为所求解，否则回到4，知道求出满足精度的点为止。

2.3.2 直接法

在无约束非线性规划中，遇到的目标函数不可导或者导函数的解析式难以表示的时候，一般使用直接法，这里介绍的是——Powell方法。

该方法主要由：基本搜索、加速搜索、调整搜索方向三部分组成

2.4 Matlab内置函数求无约束极值问题

在内置函数中求解无约束极值问题的函数有fminunc和fminsearch，用法介绍如下：

求函数的极小值

$$ min_x \quad f(x) $$

其中$x$可以是标量或者是向量。

fminunc的基本命令是：

[X,FVAL]=FMINUNC(FUN,X0,OPTIONS,P1,P2,...)

举个例子：求函数$f(x)=100(x_2-x_1^2)^2+(1-x_1)^2$的最小值、

编写m文件fun2.m，如下

function [f,g]=fun2(x);
f=100*(x(2)-x(1)^2)^2+(1-x(1))^2;  %目标函数
g=[-400*x(1)*(x(2)-x(1)^2)-2*(1-x(1));200*(x(2)-x(1)^2)];   %梯度函数

编写主函数文件并运行：

options = optimset('GradObj' , 'on' );  %优化参数
[x,y]=fminunc('fun2' ,rand(1,2),options)

即可求得函数的极小值，当然你也可以利用二阶导数：

function [f,df,d2f]=fun3(x);
f=100*(x(2)-x(1)^2)^2+(1-x(1))^2;
df=[-400*x(1)*(x(2)-x(1)^2)-2*(1-x(1));200*(x(2)-x(1)^2)];
d2f=[-400*x(2)+1200*x(1)^2+2,-400*x(1)
 -400*x(1),200];

然后主程序的优化参数记得打开海塞矩阵：

options = optimset('GradObj' , 'on' , 'Hessian' , 'on' );
[x,y]=fminunc('fun3' ,rand(1,2),options)

求多元函数的极值也可以使用fminsearch命令，其使用格式为：

[X,FVAL,EXITFLAG,IUTPUT]=FMINSEARCH(FUN,X0,OPTIONS,P1,P2,...)

例如：求函数$f(x)=sin(x)+3$取到最小值的$x$值

编写M文件fun4.m如下：

function f=fun4(x);
f=sin(x)+3;

主程序如下：

x0=2;
[x,y]=fminsearch(@fun4,x0)

即可求得在初值2附近的极小点和极小值。

3 带约束的极值问题

带有约束条件的极值问题称为约束极值问题，也叫做规划问题。

思路：将约束问题转化为无约束问题，将非线性规划转化为线性规划，主要就是把复杂变为简单。

库恩——塔克条件是非线性规划的重点，是确定某点是否为最优点的必要条件，一般来说不是充分条件（除了凸规划）

3.1 二次规划

若某非线性规划的目标函数为自变量$x$的二次函数，约束条件又全部是线性的，就称为二次规划

数学模型：

$$ min \quad \frac{1}{2}x^THx+f^Tx,\\ s.t.\quad \begin{cases} Ax \le b\\ Aeq{\cdot}x=beq\\ \end{cases} $$

这里$H$是是对称矩阵，$f,b$是列向量，$A$是相应维数的矩阵。

求解二次规划的命令是：

[X,FVAL]=QUADPROG(H,f,A,b,Aeq,beq,LB,UB,X0,OPTIONS)

返回值X是决策向量$x$的值。

举个例子：求解二次规划

$$ \begin{cases} min \quad f(x)=2x_1^2-4x_1x_2+4x_2^2-6x_1-3x_2\\ x_1+x_2 \le 3\\ 4x_1+x_2 \le 9\\ x_1,x_2 \ge 0\\ \end{cases} $$

编程如下：

h=[4,-4;-4,8];
f=[-6;-3];
a=[1,1;4,1];
b=[3;9];
[x,value]=quadprog(h,f,a,b,[],[],zeros(2,1))  %下界大于等于0，这里给了不等式约束，等式约束没有给了[]

3.2 罚函数法

利用罚函数法，可以将非线性规划的问题转化为无约束极值问题求解，这种方法被称为序列无约束最小化技术，简称SUMT。

外罚函数法：

$$ min \quad f(x)\\ s.t.\begin{cases} g_i(x) \le 0,\quad i=1,\cdots,r,\\ h_j(x) \ge 0,\quad j=1,\cdots,s,\\ k_m(x)=0,\quad m=1,\cdots,t\\ \end{cases} $$

取一个充分大的数$M \gt 0$，构造函数

$$ P(x,M)=f(x)+M\sum^r_{i=1}max(g_i(x),0)-M\sum^s_{i=1}min(h_i(x),0)+M\sum^t_{i=1}|k_i(x)| $$

以增广目标函数$P(x,M)$为目标函数的无约束极值问题

$$ min \quad P(x,M) $$

的最优解$x$也是原问题的最优解。

举个例子：求下列非线性规划

$$ \begin{cases} min \ f(x)=x_1^2+x_2^2+8\\ x_1^2-x_2 \ge 0\\ -x_1-x_2^2+2=0 \\ x_1,x_2 \ge 0\\ \end{cases} $$

解：编写M文件test.m

function g=test(x);
M=50000;   %M取到足够大
f=x(1)^2+x(2)^2+8;
g=f-M*min(x(1),0)-M*min(x(2),0)-M*min(x(1)^2-x(2),0)+...
    M*abs(-x(1)-x(2)^2+2);

在窗口运行

 [x,y]=fminunc('test',rand(2,1))

即可获得问题的解。

3.3 Matlab内置函数求解约束极值问题

这些内置函数有很多：fminbnd、fmincon、quadprog、fseminf、fminimax，之前已经介绍了函数fmincon和quadprog

3.3.1 fminbnd函数

求单变量非线性函数在区间上的极小值

$$ min_xf(x),x\in [x_1,x_2] $$

命令为：

[X,FVAL]=FMINBND(FUN,X1,X2,OPTIONS)

返回值是极小点$x$与函数的极小值，这里的FUN是m文件定义的函数或者是内置的单变量数学函数。

举个例子：求函数$f(x)=(x-3)^2-1,x\in[0,5]$的最小值

编写m文件fun5.m

function f=fun5(x);
 f=(x-3)^2-1;

在窗口输入

[x,y]=fminbnd('fun5',0,5)

即可求得极小点和极小值

3.3.2 fseminf函数

求

$$ min_x\left\{ F(x)|C(x) \le 0,Ceq(x)=0,PHI(x,w) \le 0 \right\}\\ s.t. \begin{cases} A*x \le B\\ Aeq*x=Beq\\ \end{cases} $$

其中$C(x),Ceq(x),PHI(x,w)$都是向量函数，$w$是附加的向量变量，$w$的每个分量都限定在某个区间内。

命令格式为：

X=FSEMINF(FUN,X0,NTHETA,SEMINFCON,A,B,Aeq，Beq)

举个例子：求函数$f(x)=(x_1-0.5)^2+(x_2-0.5)^2+(x_3-0.5)^2$取最小值时的$x$值，约束为：

编写M文件fun6.m定义目标函数如下：

function f=fun6(x,s);
f=sum((x-0.5).^2);

编写M文件fun7.m定义约束条件如下：

function [c,ceq,k1,k2,s]=fun7(x,s);
c=[];ceq=[];
if isnan(s(1,1))
 s=[0.2,0;0.2 0];
end
%取样值
w1=1:s(1,1):100;
w2=1:s(2,1):100;
%半无穷约束
k1=sin(w1*x(1)).*cos(w1*x(2))-1/1000*(w1-50).^2-sin(w1*x(3))-x(3)-1;
k2=sin(w2*x(2)).*cos(w2*x(1))-1/1000*(w2-50).^2-sin(w2*x(3))-x(3)-1;
%画出半无穷约束的图形
plot(w1,k1,'-',w2,k2,'+');

在命令行窗口键入：

[x,y]=fseminf(@fun6,rand(3,1),2,@fun7)  %给与初始点和约束个数

3.3.3 fminimax函数

求解

$$ min_x\left\{max_{F_i}F(x)\right\}\\ s.t.\begin{cases} A*x \le b\\ Aeq*x=Beq \\ C(x) \le 0\\ Ceq(x) = 0 \\ LB \le x \le UB\\ \end{cases} $$

其中$F(x)=\left\{F_1(x),\cdots,F_m(x)\right\}$

上述问题的命令为：

X=FMINIMAX(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON)

举个例子：求函数族$\left\{f_1(x),f_2(x),f_3(x),f_4(x),f_5(x)\right\}$取极大极小值时的$x$值。其中：

$$ \begin{cases} f_1(x)=2x_1^2+x_2^2-48x_1-40x_2+304\\ f_2(x)=-x_1^2-3x_2^2\\ f_3(x)=x_1+3x_2-18\\ f_4(x)=-x_1-x_2\\ f_5(x)=x_1+x_2-8\\ \end{cases} $$

编写M文件fun8.m定义向量函数如下：

function f=fun8(x);
f=[2*x(1)^2+x(2)^2-48*x(1)-40*x(2)+304
 -x(1)^2-3*x(2)^2
 x(1)+3*x(2)-18
 -x(1)-x(2)
 x(1)+x(2)-8];

调用函数

[x,y]=fminimax(@fun8,rand(2,1))

3.3.4 利用梯度求解约束优化问题

举个例子：已知函数$f(x)=e^{x_1}(4x_1^2+2x_2^2+4x_1x_2+2x_2+1)$,且满足非线性约束：

$$ \begin{cases} x_1x_2-x_1-x_2 \le -1.5\\ x_1x_2 \gt -10\\ \end{cases} $$

求$min_xf(x)$

题目中的目标函数梯度为：

$$ \begin{bmatrix} e^{x_1}(4x_1^2+2x_2^2+4x_1x_2+8x_1+6x_2+1)\\ e^{x_1}(4x_1+4x_2+2)\\ \end{bmatrix} $$

编写M文件fun9.m定义目标函数及梯度函数：

function [f,df]=fun9(x)
f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);
df=[exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+8*x(1)+6*x(2)+1);exp(x(1))*(4*x(2)+4*x(1)+2)];

编写M文件fun10.m定义约束条件及约束条件的梯度函数

function [c,ceq,dc,dceq]=fun10(x)
c=[x(1)*x(2)-x(1)-x(2)+1.5;-x(1)*x(2)-10];
dc=[x(2)-1,-x(2);x(1)-1,-x(1)];
ceq=[];dceq=[];

调用函数fmincon，编写主函数并运行：

%采用标准算法
options=optimset('largescale','off');
%采用梯度
options=optimset(options,'GradObj','on','GradConstr','on');
[x,y]=fmincon(@fun9,rand(2,1),[],[],[],[],[],[],@fun10,options)

非线性规划

1 非线性规划

1.1 非线性规划的Matlab解法：

1.2 非线性规划问题的迭代求解方法

2 无约束问题

2.1 一维搜索方法

2.1.1 Fibonacci法

2.1.2 0.618法

2.2 二次插值法

2.3 无约束极值问题的解法

2.3.1 解析法

2.3.1.1 梯度法（最速下降法）

2.3.1.2 Newton法

2.3.1.3 变尺度法

2.3.2 直接法

2.4 Matlab内置函数求无约束极值问题

3 带约束的极值问题

3.1 二次规划

3.2 罚函数法

3.3 Matlab内置函数求解约束极值问题

3.3.1 fminbnd函数

3.3.2 fseminf函数

3.3.3 fminimax函数

3.3.4 利用梯度求解约束优化问题

3.4 优化工具箱的用户图形界面解法