与拉格朗日乘数法比较理解



惩罚因子是一个充分大的确定的数，而拉格朗日乘子是个不确定的参数，在求最小值时确定

解最优化问题

本文版权属于重庆大学计算机学院刘骥，禁止转载



解最优化问题

最优化问题的分类
直接开始求解

downhill simplex
Conjugate Gradient
你居然不会求偏导数
Newton-CG
总结


数学原理
数学白痴的福利






最优化问题的分类

按照目标函数的可导性，可以将最优化问题分为两类： 

1. 目标函数可导，例如目标函数为$f(x)=x^2+x+1$ 

2. 目标函数不可导，例如目标函数为$f(X)=\sum_{p\in X}e^p$

对于不可导的目标函数，其解法比较特殊。采用模拟退火算法、粒子群算法、遗传算法、蚁群算法、图割（Graph Cut）算法等方法，可以求解这类问题，但通常无法得到精确的解（或者说全局最优解）。并且有些问题，仅能用特定的解法，求出特定的解。由于该类问题的复杂性，本文对不可导目标函数问题不进行讨论。如果你在实际中遇到了这类问题，可以采用如下方法解决： 

查学术论文，github上找代码；如果找不到解决方案，恭喜你，这是一个获得数学大奖的机会

本文余下部分，针对可导的目标函数进行讨论，这类问题存在大量通用的算法。人类历史发展到2017年，你甚至不需要了解具体算法的原理，也能够很happy的解决这类问题。



直接开始求解

下面我们用一个简单的最优化问题为例： 

 $$\begin{align}
\mathop{\arg\min}_{x,y}f(x,y)=x^2+y^2+1
\end{align}$$  

显然$x=0，y=0$时，$f(x,y)$有最小值1。问题在于，计算机如何求解这个问题呢？




downhill simplex

Python语言的scipy.optimize包的fmin函数实现了downhill simplex方法，也称为Nelder–Mead方法（downhill、simplex和downhill simplex是不同的概念，有兴趣可以参考维基百科）。这里不去讲解downhill simplex的原理，有需要原理的读者，请移步到https://en.wikipedia.org/wiki/Nelder–Mead_method。 

fmin函数的说明如下： 

 

我估计你看不懂，所以我不准备解释。我们直接看代码吧。

# encoding=utf-8
import scipy.optimize as opt
#定义目标函数f
def f(X):
    x=X[0]
    y=X[1]
    return x**2+y**2+1;
x=10
y=10
X=[x,y]
X=opt.fmin(f,X) #f是目标函数，X是初始猜测的解
x=X[0]
y=X[1]
print x,y 

程序执行结果如下： 

 

最终$x$和$y$的值不是0，而是一个接近0的值。这是由于数值计算毕竟是缺乏精度的。fmin函数是一个迭代算法。在执行过程中，它还同时输出了结束迭代时的最小值，迭代的次数，以及函数计算的次数。

Conjugate Gradient

downhill simplex比较简单，实现没有难度。但它的问题在于比较慢，下面我们试试 Conjugate Gradient法。使用scipy.optimize包的fmin_cg方法。这个方法需要求出目标函数的偏导数（梯度）： 

 $$\begin{align}
\frac{\partial f}{\partial x}=2x \\
\frac{\partial f}{\partial y}=2y
\end{align}$$ 




# encoding=utf-8
import scipy.optimize as opt
import numpy as np
#定义目标函数f
def f(X):
    x=X[0]
    y=X[1]
    return x**2+y**2+1;
#定义目标函数f对应的梯度函数
def gradf(X):
    x=X[0]
    y=X[1]
    gx=2*x
    gy=2*y
    return np.asarray((gx,gy)) #梯度以向量形式返回
x=10
y=10
X=[x,y]
X=opt.fmin_cg(f,X,fprime=gradf) #f是目标函数，X是初始猜测的解，增加梯度
x=X[0]
y=X[1]
print x,y 

运行结果如下： 

 

可以看出该方法比downhill simplex运行得更快。当然这是显然的，因为有了梯度信息。

你居然不会求偏导数！

OK！OK！有办法可以解决！今年是2017年，如果你还不会求偏导数，那么交给计算机吧！



X=opt.fmin_cg(f,X)

OK！计算结果如下： 

 

慢了一些，函数的演化次数增加了15次，why？因为算法用数值法来计算梯度！具体来说（你可以忽略下面的内容，把问题交给人工智能）：

def gradf(X):
    x=X[0]
    y=X[1]
    gx=(f([x+0.0000001,y])-f([x,y]))/0.0000001
    gy=(f([x,y+0.0000001])-f([x,y]))/0.0000001
    return np.asarray((gx,gy))

也就是： 

 $$\begin{align}
f'(x)=\frac{f(x+\Delta x)-f(x)}{\Delta x}
\end{align}$$ 




Newton-CG

试试Newton-CG，也就是fmin_ncg函数，该函数需要2阶偏导数构成的Hessian矩阵： 

 $$\begin{align}
\begin{bmatrix}
\frac{\partial f^2}{\partial x^2} &\frac{\partial f^2}{\partial x \partial y}\\
\frac{\partial f^2}{\partial y \partial x} &\frac{\partial f^2}{\partial y^2}
\end{bmatrix}=
\begin{bmatrix}
2&0\\
0&2
\end{bmatrix}
\end{align}$$  

代码如下：




# encoding=utf-8
import scipy.optimize as opt
import numpy as np
#定义目标函数f
def f(X):
    x=X[0]
    y=X[1]
    return x**2+y**2+1;
#定义目标函数f对应的梯度函数
def gradf(X):
    x=X[0]
    y=X[1]
    gx=2*x
    gy=2*y
    return np.asarray((gx,gy)) #梯度以向量形式返回
def hess(X):
    return np.array([[2,0],[0,2]]) #hessian矩阵
x=10
y=10
X=[x,y]
X=opt.fmin_ncg(f,X,fprime=gradf,fhess=hess) #f是目标函数，X是初始猜测的解，增加梯度
x=X[0]
y=X[1]
print x,y 

运行结果如下： 

 

非常精确！当然如果不会计算hessian，如下操作：

X=opt.fmin_ncg(f,X,fprime=gradf) #注意梯度是必须的

结果如下： 



总结

现在是2017年，基本的最优化问题已经经过了几百年的发展（不是几十年），因此有大量现成的函数库可以使用。直接调用这些库函数，你甚至不需要知道如何求导数、如何求Hessian矩阵，也可以很好的解决最优化问题（可导的）。如果不追求效率，如果精度可以接受，那么最最原始的downhill simplex方法已经能够适用于大多数场景（你需要小数点后几位的精度？）。



数学原理

大多数最优化的书籍，在讲解数学原理时，都抱着一定要让人看不懂，才能显得高深莫测的心态。其实求解可导目标函数的数学原理是极其简单的，那就是：穷举，更高效的穷举。

假设目标函数为$f(X)$，其中$X\in R^n$（注意这意味着$X$是一个向量，$f$是多元函数）。例如$X=(x,y,z)^T$，则$f(X)$表示的函数就是$f(x,y,z)$。我如此定义不是为了把你搞晕，而是为了环保！试想如果我想定义函数$f(x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8,x_9,x_{10},x_{11},x_{12},x_{13},x_{14},x_{15},x_{16},x_{17},x_{18},x_{19},x_{20},x_{21},x_{22})$，难道不是写成$f(X)$，$X\in R^{22}$更环保吗？Yes，Yes，你从来没有见过3元以上的函数，那不是因为你一直都是在做考题吗？一道考试题，要是有22个变量，你要用多大的草稿纸？现实问题则不然，考虑女性择偶问题，男方必须有钱、长得帅、有安全感、父母双亡、有车、有房、有股票、有存款、初恋、非凤凰男等等。表示为函数就是： 

 $$f(钱,颜值,安全感,父母,车,房,股票,存款,情史,是否凤凰男)$$  

学过代数吧？于是这个函数就变成了: 

 $$f(x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8,x_9,x_{10})=f(X)$$  

现在我们来讨论女性择偶函数f(X)的最优化！假定女性朋友给出的函数是这样的： 

 $$\begin{align}
f(X)=X^TC_{10\times 10}X
\end{align}$$  

这是一个二次函数（请自行展开），C是系数矩阵（我们只是假设女性朋友给出了这样的二次函数，真实情况未必）。OK，其实这个函数长什么样子，不影响后面的讨论。我们所要讨论的是$f(X)$可导的情况。如何确定$f(X)$的极值呢？


最笨的方法就是穷举！！换言之，先试试$f(X_1)$，再看看$f(X_2)$，不停的尝试$f(X_i)$，直到成为单身贵族……对，我告诉你的解法就是成为单身贵族的解法，因为解空间几乎是无限大的，穷举法太笨！

聪明一点的方式是穷举n次，找其中的最优。先试试$f(X_1)$，再看看$f(X_2)$，不停的尝试直到$f(X_n)$，找出$f(X_1)到f(X_n)$中最优的$f(X_i)$。这个方法的缺陷在于，女性朋友不一定能够找到真正的灵魂伴侣。

再聪明一点的方法固定一些变量，搜索另外一些变量。必须长得帅，像鹿晗；必须有钱，和王思聪一样；必须是初恋。你看解空间一下就变小了。接下来的搜索就简单了很多。如果找不到解，那么就修改某些固定的变量。为什么长相必须像鹿晗？郭德纲也可以接受！这样我们就可以开始新的搜索啦！事实上，这就是绝大多数人，择偶的方式。

如果变量是连续的呢？数学家发现，在这种情况下，搜索策略其实是这样的： 

1. 首先选择一个起始的解$X_0$。 

2. 比较$f(X_0+\Delta)$，如果$f(X_0+\Delta)<f(X_0)$那么$X_0+\Delta$是更好的解 

3. 另$X_1=X_0+\Delta$，跳到第1步重新执行。

上述方法的核心在于存在一个$X$的更新公式： 

 $$\begin{align}
X_{n+1}=X_{n}+\Delta_{n}
\end{align}$$ 


上述算法有三个终止条件： 

1. 迭代次数达到限额（都已经120岁的单身贵族啦！） 

2. $\Delta_{n}$趋近于0 

3. $f(X_n)-f(X_{n+1})$趋近于0

下面举个例子，目标函数如下： 

 $$\begin{align}
f(x,y)=x^2+y^2+1
\end{align}$$  

假设初始值为$X_0=(1,1)^T$，令$\Delta_n=(-0.1,-0.1)^T$，于是： 

 $$X_0=(1,1)^T\\X_1=(0.9,0.9)^T\\......\\X_{10}=(0,0)^T\\X_{11}=(-0.1,-0.1)^T$$ 


Stop！不是$X_{10}$已经达到最优了吗？为什么程序还在跑？问题在于那个$\Delta_n$，如果它能够越接近最小值越小，程序不就可以停下来了吗？越接近你的灵魂伴侣，越需要放慢搜索的速度不是吗？ 

怎么解决这个问题呢？数学家出场了。我们知道$f(X_{n+1})=f(X_n+\Delta_n)$，根据泰勒公式展开就是： 

 $$\begin{align}
f(X_{n+1})=f(X_n+\Delta_n)=f(X_n)+\frac{\partial f(X_n)}{\partial X}\Delta_n+\Delta_n^T\frac{\partial f^2(X_n)}{\partial X^2}\Delta_n /2+......
\end{align}$$  

针对该公式求$\Delta_n$的导数，并令结果等于0，则： 

 $$\begin{align}
\frac{\partial f(X_n)}{\partial X}+\frac{\partial f^2(X_n)}{\partial X^2}\Delta_n =0\\
\frac{\partial f^2(X_n)}{\partial X^2}\Delta_n=-\frac{\partial f(X_n)}{\partial X}
\end{align}$$  

现在我们可以解出$\Delta_n$！如果你对上述推导的原理感到好奇，不奇怪啊，你不是数学家对吧？总之，我们会算就可以啦。现在用新的公式来求解老问题。 

 $$f(x,y)=x^2+y^2+1\\
\frac{\partial f(x,y)}{\partial x}=2x\\
\frac{\partial f(x,y)}{\partial y}=2y\\
\frac{\partial f^2(x,y)}{\partial x^2}=2\\
\frac{\partial f^2(x,y)}{\partial x\partial y}=0\\
\frac{\partial f^2(x,y)}{\partial y^2}=2\\
\frac{\partial f(x,y)}{\partial y\partial x}=0\\
\begin{bmatrix}
2 & 0\\
0 & 2
\end{bmatrix}\Delta_n=-
\begin{bmatrix}
2x_n&2y_n
\end{bmatrix}\\
\Delta_n=-
\begin{bmatrix}
x_n&y_n
\end{bmatrix}$$  

计算过程如下： 

 $$X_0=(1,1)^T \quad \Delta_0=(-1,-1)^T\\
X_1=X_0+\Delta_0=(0,0)^T \quad \Delta_1=(0,0)^T\\
\because \Delta_1=(0,0)^T \\
\therefore X=X_1就是最优解$$ 


数学白痴表示求偏导数矩阵太难理解！！如何破！！数学家马上送上福利： 

 $$\begin{align}
\lambda\Delta_n=-\frac{\partial f(X_n)}{\partial X}
\end{align}$$  

假设$\lambda=2$，我们看看目标函数$f(x,y)=x^2+y^2+1$的计算过程：




 $$\Delta_n=-
\begin{bmatrix}
x_n&y_n
\end{bmatrix}\\X_0=(1,1)^T \quad \Delta_0=(-1,-1)^T\\
X_1=X_0+\Delta_0=(0,0)^T \quad \Delta_1=(0,0)^T\\
\because \Delta_1=(0,0)^T \\
\therefore X=X_1就是最优解$$  

运气不错！如果$\lambda=4$呢？推导如下： 

 $$\Delta_n=-
\begin{bmatrix}
0.5x_n&0.5y_n
\end{bmatrix}\\X_0=(1,1)^T \quad \Delta_0=(-0.5,-0.5)^T\\
X_1=X_0+\Delta_0=(0.5,0.5)^T \quad \Delta_1=(-0.25,-0.25)^T\\
X_2=X_1+\Delta_1=(0.25,0.25)^T \quad \Delta_2=(-0.0625,-0.0625)^T\\
......\\
若干步后收敛$$  

如果$\lambda=1$，推导如下： 

 $$\Delta_n=-
\begin{bmatrix}
2x_n&2y_n
\end{bmatrix}\\X_0=(1,1)^T \quad \Delta_0=(-2,-2)^T\\
X_1=X_0+\Delta_0=(-1,-1)^T \quad \Delta_1=(2,2)^T\\
X_2=X_1+\Delta_1=(1,1)^T \quad \Delta_2=(-2,-2)^T\\
......\\
无法收敛$$  

可见如果采用新的计算公式，$\lambda$的取值会影响算法的收敛速度，甚至无法收敛。但的确这种算法避免了计算Hessian，数学白痴非常高兴！那么它存在的问题可以解决吗？当然可以，方法就是动态调整$\lambda$的取值，从而加快收敛，避免震荡。 

最后，再来考虑一下，如果数学白痴连梯度都不会算怎么办？Nelder–Mead Method！对有高等数学背景的人而言，理解Nelder–Mead Method反而很困难！这里就不再详谈了。




数学白痴的福利

我们都不喜欢数学，还好现在是2017年，大量的工作已经由擅长数学的科学家做完了。你所需要的就是找一个算法丰富的库！比如Python的scipy或者matlab。为了写作本文，我对最优化问题进行了很多简化，首先只考虑可导的目标函数，其次也没有考虑函数变量的条件约束（比如$X$必须是整数、$X$必须大于0）。但无论问题多么复杂，这些问题已经被充分研究，充分解决。99%的情况，你只需要scipy或者matlab，剩下0.9999%的情况需要github。最后，还剩下一丢丢的未知领域，等待着学者去研究。别抢他们的饭碗，茶叶蛋已经够贵的啦！

题目：最优化问题简介

        一年多学习以来，无论是前面学习压缩感知，还是这半年学习机器学习，一直离不开最优化，比如压缩感知的基追踪类重构算法，核心问题就是一个优化问题，而机器学习中的很多算法也需要最优化的知识，比如支持向量机算法。看来必须得把最优化的基本内容学习一下了，不求理解的有多么深，至少要知道怎么用。其实前面已经写过一篇与最优化相关的内容了，就是《压缩感知中的数学知识：凸优化》这篇。
        从本篇起，开始学习一些有关最优化的基础知识，重点是了解概念和如何应用。本篇是参考文献第1.1节的一个摘编或总结，主要是把一些概念集中起来，可以随时查阅。
1、一般形式
        最优化问题的数学模型的一般形式为（以下称为最优化数学模型）：

其中，和为连续函数，通常还要求连续可微。称为决策变量，为目标函数，为约束函数，为等式约束，为不等式约束，并记等式约束的指标集为，不等式约束的指标集为。和分别是英文单词minimize(极小化)和subject
 to(受约束)的缩写。

2、概念
        如果点满足最优化数学模型中的所有约束条件就称为可行点(Feasible Point)，所有可行点的全体称为可行域(Feasible Region)，用表示。在一个可行点考虑不等式约束，如果有，就称不等式约束在点考虑不等式约束是有效约束或起作用约束(active
 constraint)，并称可行点位于约束的边界；如果有，就称不等式约束在点是无效约束或不起作用约束(inactive
 constraint)；对于一个可行点，如果没有一个不等式约束是有效的，就称是可行域的内点，不是内点的可行点就是可行域的边界点。显然在边界点至少有一个不等式约束是有效约束，当存在等式约束时，任何可行点都要满足等式约束，因此不可能是等式约束的内点。

        如果一个可行点满足，则称为最优化问题的全局最优解（或总体最优解）；如果可行点满足，则称为最优化问题的严格全局最优解（或严格总体最优解）；对于可行点，如果存在一个邻域使得成立，则称为最优化问题的局部最优解，其中是一个小的正数；对于可行点，如果存在一个邻域使得成立，则称为最优化问题的严格局部最优解。如下图所示，点是严格全部极小解，和则是局部极小解，其中是严格局部极小解，而是非严格局部极小解。（注：附近有一段线是水平的）



        一般常见的最优化方法只适用于确定最优化问题的局部最优解，有关确定全局最优解的最优化方法属于最优化问题的另一个领域——全局最优化。然而，如果最优化问题的目标函数是凸的，而可行域是凸集，则问题的任何最优解（不一定唯一）必是全局最优解，这样的最优化问题又称为凸规划。进一步，对于凸集上的严格凸函数的极小化问题，存在唯一的全局最优解。
3、分类
1）约束最优化问题
        只要在问题中存在任何约束条件，就称为约束最优化问题。
        只有等式约束时，称为等式约束最优化问题，数学模型为：

        只有不等式约束时，称为不等式约束最优化问题，数学模型为：




        既有等式约束，又有不等式约束，则称为混合约束优化问题（或一般约束优化问题）；
        把简单界约束优化问题称为盒式约束优化问题（或有界约束优化问题），数学模型为：


2）无约束最优化问题
        如果问题中无任何约束条件，则称为无约束最优化问题，数学模型为：


3）连续与离散最优化问题
        决策变量的取值是连续的，称为连续最优化问题；
        决策变量的取值是离散的，称为离散最优化问题，又称为组合最优化问题。如整数规划、资源配置、邮路问题、生产安排等问题都是离散最优化问题的典型例子，求解难度比连续最优化问题更大。
4）光滑与非光滑最优化问题
        如果最优化数学模型中的所有函数都连续可微，则称为光滑最优化问题；
        只要有一个函数非光滑，则相应的优化问题就是非光滑最优化问题。
5）线性规划问题
        对于连续光滑最优化问题，如果最优化数学模型中的所有函数都是决策变量的线性函数，则称为线性规划问题。线性规划问题的一般形式为：


其中。矩阵向量形式为



其中，
，


6）二次规划问题
        对于连续光滑最优化问题，如果最优化数学模型中的目标函数是决策变量的二次函数，而所有约束都是决策变量的线性函数，则称为二次规划问题。二次规划问题的一般形式为：

其中，为纯量，为阶对称矩阵。如果为半正定矩阵，则称此规划为凸二次规划，否则为非凸规划。对于非凸规划，由于存在比较多的驻点，求解比较困难。
7）非线性最优化问题
        只要最优化数学模型中的函数有一个关于决策变量是非线性的，则称为非线性最优化问题。
        非线性最优化问题是最一般的最优化问题，而线性规划和二次规划问题却是相当重要的特殊的最优化问题，因为在实际中形成的许多最优化问题都是线性规划问题或二次规划问题，而且在用迭代法求非线性最优化问题的最优解时我们常常用线性规划或二次规划来局部近似原非线性最优化问题，并通过求所得近似问题的最优解来对已有最优解的估计进行改进。
 
参考文献：
孙文瑜, 徐成贤, 朱德通.最优化方法(第二版)[M]. 北京:高等教育出版社, 2010.

MATLAB 求解最优化问题



MATLAB 优化工具箱解线性规划



模型1

minz=cXs.t.AX≤b
命令：x=linprog(c,A,b)



模型2

minz=cXs.t.AX≤bAeq⋅X=beq
命令：x=linprog(c,A,b,Aeq,beq)

注意：若没有不等式：AX≤b存在，则令A=[ ]，b=[ ]



模型3

minz=cXs.t.AX≤bAeq⋅X=beqVLB≤X≤VUB
命令：[1] x=linprog(c,A,b,Aeq,beq,VLB,VUB)，[2] x=linprog(c,A,b,Aeq,beq,VLB,VUB,x0)

注意：[1] 若没有等式约束：Aeq⋅X=beq，则令Aeq=[ ]，beq=[ ]；[2] 其中x0表示初始点



命令：[x,fval]=linprog(…)

返回最优解x及x出的目标函数的值fval



求解优化问题的主要函数



1. MATLAB求解优化问题的主要函数

类型
  
模型
  
基本函数名
一元函数极小
  
min F(x)s.t. x1<x<x2
  
x=fminbnd(′F′,x1,x2)
无约束极小
  
min F(X)
  
X=fminunc(′F′,x0)X=fminsearch(′F′,x0)
线性规划
  
min cTXs.t. AX≤b
  
X=linprog(c,A,b)
二次规划
  
min 12XTHX+cTXs.t. AX≤b
  
X=quadprog(H,c,A,b)
约束极小（非线性规划）
  
min F(X)s.t. G(X)≤b
  
X=fmincon(′FG′,x0)
达到目标问题
  
min rs.t. F(x)−wr≤goal
  
X=fgoalattain(′F′,x,goal,w)
极小极大问题
  
min maxx {Fi(x)}s.t. G(x)≤0
  
X=fminimax(′FG′,x0)
2. 优化函数的输入变量

变量
  
描述
  
调用函数
f
  
线性规划的目标函数f∗X或二次规划的目标函数X′∗H∗X+f∗X中线性项的系数向量
  
linprog, quadprog
fun
  
非线性优化的目标函数。fun必须为行命令对象或M文件、嵌入函数或MEX文件的名称
  
fminbnd, fminsearch,fminunc, fmincon, lsqcurvefit, lsqnonlin, fgoalattain, fminimax
H
  
二次规划的目标函数X′∗H∗X+f∗X中二次项的系数矩阵
  
quadprog
A, b
  
A矩阵和b向量分别为线性不等式约束：AX≤b中的系数矩阵和右端向量
  
linprog, quadprog, fgoalattain, fmincon, fminimax
Aeq, beq
  
Aeq矩阵和beq向量分别为线性等式约束：Aeq⋅X=beq中的系数矩阵和右端向量
  
linprog, quadprog, fgoalattain, fmincon, fminimax
vlb, vub
  
X的下限和上限向量：vlb≤X≤vub
  
linprog, quadprog, fgoalattain, fmincon, fminimax, lsqcurvefit, lsqnonlin
X_0
  
迭代初始点坐标
  
除fminbnd外所有优化函数
X1, X2
  
函数最小化的区间
  
fminbnd
options
  
优化选项参数结构，定义用于优化函数的参数
  
所有优化函数
3. 优化函数的输出变量表

变量
  
描述
  
调用函数
x
  
由优化函数求得的值。若exitflag>0，则x为解；否则x不是最终解，它只是迭代终止时优化过程的值
  
所有优化函数
fval
  
解x处的目标函数值
  
linprog, quadprog, fgoalattain, fmincon, fminimax, lsqcurvefit, lsqnonlin, fminbnd
exitflag
  
描述退出条件：exitflag>0，表明目标函数收敛于解x处；exitflag0=，表明目标函数评价或迭代的最大次数；exitflag<0，表明目标函数不收敛
  

output
  
包含优化结果信息的输出结构：Iterations：迭代次数；Algorithm：所采用的算法；FuncCount：函数评价次数
  
所有优化函数
4. 控制参数options的设置



用MATLAB解无约束问题



1. 一元函数无约束优化问题

min f(x)s.t. x1≤x≤x2
常用格式如下



x=fminbnd(fun,x1,x2);
x=fminbnd(fun,x1,x2,options);
[x,fval]=fminbnd(...);
[x,fval,exitflag]=fminbnd(...);
[x,fval,exitflag,output]=fminbnd(...);
%%其中3、4、5的等式右边可用1或2的等式右边

函数fminbnd的算法基于黄金分割法和二次插值法，它要求目标函数必须是连续函数，并可能只给出局部最优解。



2. 多元函数无约束优化问题

标准型为：min F(X)

命令格式为：



x=fminunc(fun,X0);  %或x=fminsearch(fun,X0)
x=fminunc(fun,X0,options);  %或x=fminsearch(fun,X0,options)
[x,fval]=fminunc(...);  %或[x,fval]=fminsearch(...)
[x,fval,exitflag]=fminunc(...); %或[x,fval,exitflag]=fminsearch(...)
[x,fval,exitflag,output]=fminunc(...);  %或[x,fval,exitflag,output]=fminsearch(...)

fminsearch是用单纯性法寻优

fminunc的算法：

fminunc为无约束优化提供了大型优化和中型优化算法。 
由options中的参数LargeScale控制：LargeScale='on'使用大型算法，LargeScale='off'使用小型算法
fminunc为中型优化算法的搜素方向提供了4种算法，由options中的参数HessUpdate控制 
HessUpdate='bfgs'（默认值），拟牛顿法的BFGS公式
HessUpdate='dfp'，拟牛顿法的DFP公式
HessUpdate='steepdesc'，最速下降法
fminunc为中型优化算法的步长一维搜索提供了两种算法，由options中参数LineSearchType控制： 
LineSearchType='quadcubic'（缺省值），混合的二次和三次多项式插值
LineSearchType='cubicpoly'，三次多项式插值
使用fminunc和fminsearch可能会得到局部最优解



用MATLAB解非线性规划

标准型为： 
min F(X)s.t. AX≤bAeq⋅X=beqG(X)≤0Ceq(X)=0VLB≤X≤VUB


其中X为n维变元向量，G(X)与Ceq(X)均为非线性函数组成的向量，其他变量的含义和线性规划、二次规划相同。MATLAB求解：

首先建立M文件fun.m，定义目标函数F(X)

function f=fun(X);
f=F(X);
若约束条件中有非线性约束：G(X)≤0或Ceq(X)=0，则建立M文件nonlcon.m定义函数G(X)与Ceq(X)：

function [G,Ceq]=nonlcon(X);
G=...
Ceq=...
建立主程序，非线性规划求解的函数是fmincon，命令的基本格式如下：

x=fmincom('fun',X0,A,b);
x=fmincon('fun',X0,A,b,Aeq,beq);
x=fmincon('fun',X0,A,b,Aeq,beq,VLB,VUB);
x=fmincon('fun',X0,A,b,Aeq,beq,VLB,VUB,'nonlcon');
x=fmincon('fun',X0,A,b,Aeq,beq,VLB,VUB,'nonlcon',options);
[x,fval]=fmincon(...);
[x,fval,exitflag]=fmincon(...);
[x,fval,exitflag,output]=fmincon(...);
fmincon函数提供了大型优化算法和中型优化算法。默认时，若在fun函数中提供了梯度（options参数的GradObj设置为on），并且只有上下界存在或只有等式约束，fmincon函数将选择大型算法。当既有等式约束又有梯度约束时，使用中型算法。

fmincon函数的中型算法使用的是序列二次规划法。在每一步迭代中求解二次规划子问题，并用BFGS法更新拉格朗日Hessian矩阵。

fmincon函数可能会给出局部最优解，这与初值X0的选取有关。