数学建模（三）：非线性规划问题

&1.非线性规划
1.当一个规划问题的约束条件或者是目标函数中包含至少一个非线性函数，那么这种规划问题就是⎡非线性规划问题⎦。
举个例子说明：
设决策变量x_i的意义为对第i个项目投资与否，则x_i=1代表对这个项目投资，x_i=0代表对这个项目不投资；可以得到对于投资金额来说，由于至少要对一个项目投资，并且总投资金额不能超过总资金A

对于决策变量，它是一个只能二元取值的变量，所以有x_i(1-x_i)=0,i=1,2,…,n
我们这里对于投资收益最大化的理解是投资总收益与总投资的比值最大化，所以这个规划问题综上所述可以表示为

《数学建模算法与应用》中还给出了对于一个实际问题把它转化为非线性规划问题的注意要点，但过于泛泛而谈，以至于光看这些注意要点对于所有的规划问题都能成立。我认为非线性规划问题在建模时应当注意⎡是否有必要引入非线性的目标函数和约束条件⎦，一个线性规划问题总是比非线性规划问题容易求得最优解，因为线性规划问题的最优解总是会在问题的可行域的边界上，而非线性规划问题的最优解却可以存在于问题的可行域中的任意一点；
非线性规划是可以通过matlab的FMINCON函数求解的；

2.如何求解非线性规划问题？
一种方法是迭代法。在讲述迭代法之前，我们需要知道整体最优解、严格整体最优解、局部最优解和严格局部最优解的定义，这里假定你已经知道了这些概念代表的含义。使用迭代法的基本思想是：选定初始点x⁰∈Rⁿ，按照某种特定的迭代规则，可以由x⁰产生一系列的解，由它们组成点集{x^k}：
1）当{x^k}中元素个数是有限时，那该点集的最优一个点就是该非线性规划问题的最优解；
2）当{x^k}中元素个数是无限时，那么它有极限点，并且这个极限点就是非线性规划问题的最优解；
假设x^k∈Rⁿ是第k轮迭代点，那么在它不是NP问题（Non-linear Programming，非线性规划）最优解情况下，它必定有下一个迭代点x^k+1，它们之间的关系满足：

x^k+1=x^k+t_kp^k

式中的t_k是一个实数，p^k∈Rⁿ，||p^k||=1，并且p^k方向为从x^k向x^k+1的方向。这里的p^k我们称它为第k轮搜索方向，t_k我们称它为沿p^k方向的步长。
我们引入两个概念下降方向和可行方向：
1）下降方向：指当存在大于0的步长，使得目标函数在第k+1轮迭代中的结果小于第k轮迭代结果，此时的p^k方向称为目标函数在x^k点的下降方向；
2）可行方向：指当存在大于0的步长，使得第k+1轮的迭代结果在NP问题的可行域K中，此时的p^k方向称为x^k点关于可行域K的可行方向；
既是可行方向又是下降方向，我们称之为可行下降方向；接下来给出基本迭代格式求解NP问题的一般步骤：
0^。选取初始点x⁰，令k=0；
1^。构造搜索方向，按照迭代规则，得到可行下降方向p^k；
2^。寻找合适的步长t_k，求出x^k+1;若满足终止条件，则停止迭代；
3^。以x^k+1代替x^k，返回1^。步；

3.凸函数、凸规划
凸函数、严格凸函数的定义我们在此略去，考虑非线性规划
假定f(x)为凸函数，g_j(x)为凸函数，那么这样的规划称为凸规划。凸规划的可行域一定是凸集，它的局部最优解就是全剧最优解，它的最优解的集合可以形成一个凸集。当目标函数f(x)是严格凸函数时，假如凸规划问题有最优解，它的最优解是唯一的。

&2.无约束问题
1.一维搜索方法
沿某一已知方向求目标函数的极小值。常用的一维搜索方法：（1）试探发（Fibonacci法，黄金比例法）；（2）插值法（抛物线插值法，三次插值法）；（3）微积分中的求根法（切线法，二分法等）
在讲述详细的搜索方法之前，我们先考虑一维极小化问题，假设函数f(t)在a≤t≤b范围内有唯一的极值，那么我们通过不断缩短区间，来寻求近似最优解t^*；在[a,b]区间中任取两个关于(a+b)/2对称的点t₁,t₂，假定t₁<t₂，通过测定法f(t₁)和f(t₂)的值，我们可以得到：
1）若f(t₁)<f(t₂)，则必有t^*∈[a,t₂]；
2）若f(t₂)<f(t₁)，则必有t^*∈[t₁,b]；
3）若f(t₁)=f(t₂)，则必有t^*∈[t₁,t₂]；
(4.22 后续待更新)
1.1Fibonacci法
引入Fibonacci数列定义，假如数列{F_n}满足关系

F₀=F₁=1
F_n=F_n-2+F_n-1,n=2,3,…,

数列{F_n}就被称为Fibonacci数列，F_n称为第n个Fibonacci数，相邻两个Fibonacci数的比值F_n-1/F_n称为Fibonacci分数，Fibonacci分数始终是小于等于1的。
接下来我们介绍通过不断缩小区间[a,b]，求得单峰目标函数f(t)在这个区间的极值的Fibonacci法：
1^。给定搜索精度𝛅>0，通过关系式(b-a)/F_n≤𝛅，可以得到n的值；
2^。k=1，根据t₁=a+F_n-2/F_n*(b-a),t₂=a+F_n-1/F_n (b-a)；可以计算出最初两个搜索点；
3^。当k<n-1时，计算函数值f(t₁)与f(t₂)，若f(t₁)<f(t₂)，则令b=t₂;t₂=t₁;t₁=b+F_n-1-k/F_n-k(a-b)，即左端点a不动，更新其他三个参数；
若f(t₁)>f(t₂)，则令a=t₁;t₁=t₂;t₂=a+F_n-1-k/F_n-k *(b-a)，即右端点b不动，更新其他三个参数；
k += 1;
4^。重复3^。，直到k=n-1；
5^。此时t₁=t₂=(b+a)/2，此时令t₁=(b+a)/2，t₂=a+(1/2+𝛆)(b-a)；𝛆为任意小的数；取此时f(t₁)与f(t₂)中的较小值作为近似极小值，相应的区间为[a,t₂]或[t₁,b]

1.2黄金比例法
一种Fibonacci法的近似方法，实现起来较为容易，黄金比例指一个事物的(√5-1)/2倍，约为0.618，令𝛚=(√5-1)/2，这个数满足𝛚/1=(1-𝛚)/𝛚，而且它是Fibonacci分数当n趋近于无穷时的极限；它和Fibonacci方法较为相似，不过缩短区间的比例为0.618，也就是𝛚的近似数；

2.二次插值法
对于f(t)在[a,b]上连续时的极小化问题，可以考虑使用多项式插值进行求解。通过不断用低次多项式（指二次以及二次以下）来近似f(t)，用插值多项式的极小点来逼近原问题的最优解；

3.无约束极值问题
当不限制区间，也不限制极值个数时，这样的极值问题我们称它为无约束极值问题，这种问题有解析法和直接法两种方法，区别是前者会用到导数，而后者不需要；我们分别对这两种方法进行介绍

3.1解析法
3.1.1梯度法（最速下降法）
我们在之前迭代法中介绍了最基本的迭代格式x^k+1=x^k+t_kp^k，我们介绍p^k时，说明了它是一个方向，能够使得函数f(x)的值在点x^k下降；假设我们要求这个方向p^k上f(x)的值在点x^k下降得最快，那么由微积分知识告诉我们这个方向应该是f(x)在点x^k的负梯度方向（关于这一点的说明可以参考这里，我并不打算花很多篇幅来说明此事，但我希望本文读者假如无法解决⎡为什么函数在某一点的最速下降方向是它的负梯度方向？⎦这样的问题，应该去仔细读一读这篇文章），所以我们令p^k=-∇f(x^k)；梯度法的终止条件是∇f(x^x)=0或||∇f(x^k)||≤𝛆，我们来谈一谈具体步骤：
1^。选取初始点x⁰，给定终止误差𝛆，令k=0；
2^。求梯度向量，计算∇f(x^k)，当∇f(x^x)=0或||∇f(x^k)||≤𝛆时，停止迭代，输出x^k。
3^。构造p^k=-∇f(x^k)，并求得t_k，使得f(x^k+t_kp^k) = min f(x^k+tp^k)（t≥0）
4^。令x^k+1=x^k+t_kp^k，k+=1，转2^。
关于如何求梯度，可以自行解决，不做赘述。

3.1.2Newton法（在讲述这一块时，我假定你已经懂得了多元函数的泰勒展开以及黑塞矩阵及黑塞矩阵正定的性质）
对于一个多元目标函数f(x)而言，在点x^k处通过泰勒展开可以获得它的二次逼近式：
f(x)≈Q(x)=f(x^k)+∇f(x^k)^T(x-x^k)+1/2*(x-x^k)^T∇²f(x^k)(x-x^k)；我们假设黑塞矩阵∇²f(x^k)是正定的，因此Q(x)的驻点即为它的极值点。若设x^k+1为Q(x)的极值点，有∇Q(x^k+1)=∇f(x^k)+∇²f(x^k)(x^k+1-x^k)=0，就可以得到x^k+1=x^k-[∇²f(x^k)]^-1∇f(x^k)，此时[∇²f(x^k)]^-1∇f(x^k)被我们称作Newton方向，并且此时步长t_k取固定值1，接下来我们给出Newton法具体步骤：
1^。选取初始点x⁰，给定终止误差𝛆>0，令k=0；
2^。求梯度向量，即∇f(x^k)，若||∇f(x^k)||≤𝛆，则停止迭代，输出x^k。否则进行3^。；
3^。计算[∇²f(x^k)]^-1，取p^k=-[∇²f(x^k)]^-1∇f(x^k)；
4^。求得下一迭代点x^k+1=x^k+p^k；
Newton法的优点是收敛速度快；缺点是初始点对于求解过程的影响较大，选择了一个不大理想的初始点，可能会导致求不出解。

3.1.3 变尺度法（Variable Metric Algorithm）
变尺度法避免了计算二阶导数矩阵及其求逆过程，又比梯度法收敛速度快，使得它获得了很高的声誉。我们接下来介绍它一种方法——DFP法。
为了不计算Newton方法中的二阶导数矩阵及其逆矩阵，我们使用一个构造矩阵，它能够近似于[∇²f(x^k)]^-1，因此这种方法也叫做拟牛顿法；
关于矩阵构造过程，我们不多赘述，我们要清楚的是这是一个迭代过程，每一步中得到的近似矩阵为H^(k)，而迭代关系满足H^(k+1)=H^(k)+ΔH^(k)，ΔH^(k)被叫做校正矩阵；经过一些数学过程，我们可以得到

这里的ΔG^(k)=∇f(x^k+1)-∇f(x^k)；
于是我们得到确切的迭代关系，而一般我们取H⁽⁰⁾为单位矩阵，下面我们给出DFP法的步骤：
1^。给定初始点x⁰以及终止误差𝛆；
2^。若||∇f(x⁰)||≤𝛆，则停止迭代，输出x^k。否则进行3^。；
3^。令H⁽⁰⁾=I(单位矩阵)，p⁰=-H⁰∇f(x⁰)，在p⁰方向进行一维搜索，确定最佳步长t₀使得f(x⁰+t₀p⁰)最小，并得到x¹=x⁰+t₀p⁰;
4^。计算[∇²f(x^k)]^-1，||∇f(x^k)||≤𝛆，则停止迭代，输出x^k，否则计算H^(k)，并令p^k = -H^(k)∇f(x^k)，得到x^k+1=x^k+t_kp^k；
5^。若x^k+1满足精度要求，则得到近似解x^k+1，否则转会4^。，直到满足要求；

3.2 直接法
在函数不可导或者是导函数的解析式难以表示时，人们一般需要使用直接搜索方法。我们介绍其中一种Powell方法。
这个方法主要由基本搜索、加速搜索和调整搜索方向三部分组成，具体步骤如下：
1^。选取初始点x⁰以及n个线性无关初始方向，组成P={p⁰,p¹,p²,…,p^n-1}，给定终止误差𝛆>0，令k=0；
2^。进行基本搜索，令y⁰=x^k，沿着P中的方向进行一维搜索，对应的可以得到y¹,y²,…,yⁿ，即y^j=y^j-1+t_j-1p^j-1;f(y^j)=min f(y^j-1+tp^j-1)；
3^。构造加速方向。令pⁿ=yⁿ-y⁰，若||pⁿ||≤𝛆，停止迭代，输出x^k+1=yⁿ，否则进行4^。
4^。确定调整方向，有f(y^m-1)-f(y^m)=max{f(y^j-1)-f(y^j)|1≤j≤n}找到m,当f(y⁰)-2f(yⁿ)+f(2yⁿ-y⁰)<2[f(y^m-1)-f(y^m)]成立进行5^。，都则进行6^。。
5^。调整搜索方向组，x^k+1=yⁿ+t_npⁿ：f(yⁿ+t_npⁿ)=min f(yⁿ+tpⁿ)
同时，更新P为{p⁰,p¹,…p^m-1,p^m+1,…,p^n-1,pⁿ}，k+=1,转2^。；
6^。不调整搜索方向组，令x^k+1=yⁿ,k+=1，转2^。；

&3 约束极值问题
求解约束极值问题会比求解无约束极值问题要难得多，可以采用的方法有：将约束问题转化为无约束问题；将非线性规划问题化为线性规划问题，以及能将复杂问题转化为简单问题的其他方法；
1）罚函数法
利用该方法可以讲非线性规划问题求解转化为求解一系列无约束极值问题，因而这种方法又称为序列无约束最小化技术（SUMT,Sequential Unconstrained Minization Technique）。
罚函数法求解的思想是，利用问题中的约束函数做出适当的罚函数，由此构造出带参数的增广目标函数，把问题转化为无约束非线性规划；罚函数有外罚函数法和内罚函数法，下面介绍外罚函数法；
举个例子

我们选取一个充分大的正数M，构造函数

这里G(x)=[g₁(x) … g_r(x)]，H(x)=[h₁(x) … h_s(x)],K(x)=[k₁(x) … k_t(x)],则我们可以得到一个无约束极值问题min P(x,M)。