非线性规划

作者: Acapella_Zhang | 来源:发表于2019-05-01 12:17 被阅读0次

1. 非线性规划

1.1 示例以及定义

如果目标函数或约束条件中包含非线性函数，就称这种规划问题为非线性规划问题。一般说来，解非线性规划要比解线性规划问题困难得多。而且，也不象线性规划有单纯形法这一通用方法，非线性规划目前还没有适于各种问题的一般算法，各个方法都有自己特定的适用范围。

1.2 线性规划与非线性规划的区别

如果线性规划的优解存在，其优解只能在其可行域的边界上达到（特别是可行域的顶点上达到）；而非线性规划的优解（如果优解存在）则可能在其可行域的任意一点达到。

1.3非线性规划的Matlab 解法

Matlab 中非线性规划的数学模型写成以下形式

$minf(x) \left\{\begin{array}{l}{A x \leq B} \\ {A e q \cdot x=B e q} \\ {C(x) \leq 0} \\ {Ceq(x)=0}\end{array}\right.$
其中 $C(x)和Ceq(x)$ 是非线性函数
Matlab 中的命令是
X=FMINCON(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS)
NONLCON 是用 M 文件定义的非线性向量函数 $C(x)和Ceq(x)$
给出例子如下
$\min f(x)=x_{1}^{2}+x_{2}^{2}+x_{3}^{2}+8$
$\begin{array}{l}{x_{1}^{2}-x_{2}+x_{3}^{2} \geq 0} \\ {x_{1}+x_{2}^{2}+x_{3}^{3} \leq 20} \\ {-x_{1}-x_{2}^{2}+2=0}\\{x_{2}+2 x_{3}^{2}=3\\ x_{1}, x_{2}, x_{3} \geq 0}\end{array}$

function f=fun1(x); %定义目标函数
f=sum(x.^2)+8;

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]; %非线性等式约束

求解程序：

[x,y]=fmincon('fun1',rand(3,1),[],[],[],[],zeros(3,1),[],'fun2')
>>
x =

    0.5522
    1.2033
    0.9478


y =

   10.6511

2.无约束问题

2.1 一维搜索方法

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

试探法（“成功—失败”，斐波那契法，0.618 法等）
插值法（抛物线插值法，三次插值法等）
微积分中的求根法（切线法，二分法等）

下面有两个通过不断地缩短区间[a,b]的长度，来搜索得到近似优解的两个方法。

2.1.1 斐波那契法

2.1.2 0.618法

2.2 无约束极值问题的解法

迭代法大体上分为两点：一是用到函数的一阶导数或二阶导数，称为解析法。另一是仅用到函数值，称为直接法。

2.3.1 解析法

梯度法

其中 $t^k$ 为步长，通过求在这点去最小函数值时的 $t^k$
给出一个例子如下：
【例】
用最速下降法求解 $min f(x) = x_1^2+25x_2^2$
要求初始点为 $x^0 = [2,2]^T$

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);
while norm(g)>0.000001
    p = -g/norm(g);
    t = 1.0;
    f = detaf(x+t*p);
    while f>f0   %寻找让函数值最小的t
        t = t/2;
        f = detaf(x+t*p);
    end
    x = x+t*p;
    [f0,g] = detaf(x);
end
x,f0

牛顿法
考虑目标函数在点 $x^k$ 的二次逼近
$f(x) \approx Q(x)=f\left(x^{k}\right)+\nabla f\left(x^{k}\right)^{T}\left(x-x^{k}\right)+\frac{1}{2}\left(x-x^{k}\right)^{T} \nabla^{2} f\left(x^{k}\right)\left(x-x^{k}\right)$
假定海塞矩阵正定：
$\nabla^{2} f\left(x^{k}\right)=\left[ \begin{array}{ccc}{\frac{\partial^{2} f\left(x^{k}\right)}{\partial x_{1}^{2}}} & {\cdots} & {\frac{\partial^{2} f\left(x^{k}\right)}{\partial x_{1} \partial x_{n}}} \\ {\vdots} & {\ldots} & {\vdots} \\ {\frac{\partial^{2} f\left(x^{k}\right)}{\partial x_{n} \partial x_{1}}} & {\cdots} & {\frac{\partial^{2} f\left(x^{k}\right)}{\partial x_{n}^{2}}}\end{array}\right]$

2.4 matlab求解无约束极值

3. 约束极值问题

带有约束条件的极值问题称为约束极值问题，也叫规划问题。求解约束极值问题要比求解无约束极值问题困难得多。为了简化其优化工作，可采用以下方法：将约束问题化为无约束问题；将非线性规划问题化为线性规划问题，以及能将复杂问题变换为较简单问题的其它方法。库恩—塔克条件是非线性规划领域中重要的理论成果之一，是确定某点为优点的必要条件，但一般说它并不是充分条件（对于凸规划，它既是优点存在的必要条件，同时也是充分条件）。

3.1 二次规划

若某非线性规划的目标函数为自变量 x的二次函数，约束条件又全是线性的，就称这种规划为二次规划。
$\min \frac{1}{2} x^{T} H x+f^{T} x$
$s.t.\left\{\begin{array}{l}{A x \leq b} \\ {A e q \cdot x=b e q}\end{array}\right.$
【例】求解如下的例子

【注意】要提出 $\frac{1}{2}$

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)) 
>>
x =

    1.9500
    1.0500


value =

  -11.0250

3.2 罚函数法

利用罚函数法，可将非线性规划问题的求解，转化为求解一系列无约束极值问题，因而也称这种方法为序列无约束小化技术
罚函数法求解非线性规划问题的思想是，利用问题中的约束函数作出适当的罚函数，由此构造出带参数的增广目标函数，把问题转化为无约束非线性规划问题。主要有两种形式，一种叫外罚函数法，另一种叫内罚函数法，下面介绍外罚函数法。
$min f(x) \left\{\begin{array}{l}{g_{i}(x) \leq 0, i=1, \cdots, r} \\ {h_{j}(x) \geq 0, j=1, \cdots, \mathrm{s}} \\ {k_{m}(x)=0, m=1, \cdots, t}\end{array}\right.$
取一个充分大的数M>0，构造函数如下：
$P(x, M)=f(x)+M \sum_{i=1}^{r} \max \left(g_{i}(x), 0\right)-M \sum_{i=1}^{s} \min \left(h_{i}(x), 0\right)+M \sum_{i=1}^{T}\left|k_{i}(x)\right|$
直观上可以理解为若 $g(x)>0$ 则给这个目标函数一个很大的惩罚，而如果 $g(x)>0$ 则对该目标函数无影响.
或者写成：

在matlab中可以直接使用min，max，sum函数。问题转化为增广目标函数的无约束优化问题，然后用fminunc()函数解决

【例】
$\left\{\begin{array}{c}{\min f(x)=x_{1}^{2}+x_{2}^{2}+8} \\ {x_{1}^{2}-x_{2} \geq 0} \\ {-x_{1}-x_{2}^{2}+2=0} \\ {x_{1}, x_{2} \geq 0}\end{array}\right.$

function g=test1(x);
M=50000;
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);

或者写成矩阵形式：

function g=test2(x);
M=50000;
f=x(1)^2+x(2)^2+8;
g=f-M*sum(min([x';zeros(1,2)]))-M*min(x(1)^2-x(2),0)+M*abs(-x(1)-x(2)^2+2);

其中的min([x';zeros(1,2)])表示 $x_1,x_2$ 都大于0
再在命令行中输入

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

    1.3118
    0.8296


y =

   10.4090

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

    1.1810
    0.9050


y =

   10.2140

可以看出两次的效果略有偏差。但是我们不满足于此，由于问题的规模较小，尝试使用fmincon求得精确得最优解

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

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

[x,y] = fmincon('fun1',rand(2,1),[],[],[],[],zeros(3,1),[],'fun2',options)
>>
x =

    0.5000
    1.2247


y =

    9.7500

可以看出罚函数法虽然收敛速度快，但是精确度不是很高。当问题规模较大，不好求解时，可以考虑使用。

3.3 使用梯度法求解约束线性规划

$f(x)=e^{x_{1}}\left(4 x_{1}^{2}+2 x_{2}^{2}+4 x_{1} x_{2}+2 x_{2}+1\right)$
其中约束条件为：
$\left\{\begin{array}{l}{x_{1} x_{2}-x_{1}-x_{2} \leq-1.5} \\ {x_{1} x_{2} \geq-10}\end{array}\right.$
当使用梯度求解上述问题时，效率更高并且结果更准确。题目中目标函数的梯度为（对 $x_1,x_2$ 分别求偏导数）：
$\left[ \begin{array}{l}{e^{x_{1}}\left(4 x_{1}^{2}+2 x_{2}^{2}+4 x_{1} x_{2}+8 x_{1}+6 x_{2}+1\right)} \\ {e^{x_{1}}\left(4 x_{1}+4 x_{2}+2\right)}\end{array}\right]$

编写fun10定义目标函数和梯度函数：

function [f,df]=fun10(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)];

编写fun11定义约束条件，以及约束条件的梯度函数

function [c,ceq,dc,dceq]=fun11(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=[];

调用fincon()

options=optimset('GradObj','on','GradConstr','on');%采用梯度
[x,y]=fmincon(@fun10,rand(2,1),[],[],[],[],[],[],@fun11,options)
>>
x =

   -9.5473
    1.0474


y =

    0.0236

3.5 optimtool工具箱的使用

在命令行中输入optimtool可以打开图形界面

对于上一个问题，只要选择方法（solver）后，在相应位置输入@fun10，@fun11，点击start就可以