《数学建模方法与分析》（《Mathematical Modeling》）——最优化模型

这学期要学最优化方法，以及要搞数学建模，所以这本圣经必须得安排了。

本文不能替代书本，只是我个人的结构化总结，还会加上个人的理解，结合不同领域的知识，所以谨慎看。

我个人感觉，这本书应该放在大一下，和线代数分下同步进行，一边学，一边练，正好，再不济也得大二上。但是我大二下才学，要是早点看透学校的本质，早点自学就好了。

对了，另一本叫《数学建模》的翻译本，例子太多了，不够精悍，相比来说，这本书的例子更多地集中于课后题，而章节里都似乎例题，比较灵活。

最优化模型

单变量最优化

常见的问题有三要素

1. 决策变量
2. 约束条件
3. 目标函数

我们要在约束条件的限制下，找出最优的决策变量组合，来使得目标函数达到最大或者最小。

五步方法

五步法大致是先将问题转化为数学模型，然后求解数学模型，最后再将数学答案转换成自然答案。

1. 提出问题
2. 选择建模方法
3. 推导模型的数学表达式
4. 求解模型
5. 回答问题

提出问题

• 首先列出所有变量

变量的名字应该具有实际意义，比如时间t（time），重量w（weight），成本c（cost），毛利润r（reward），净利润p（pure）等等。

同时注意不要混淆变量和常量，常量可以使用大写来区分。

变量应该具有明确的单位，否则没有意义。

• 然后列出所有假设

就是我们的约束条件，以及各种已知的东西。因为前面已经列出所有变量，所以假设就用前面的变量去描述

• 最后写出目标

请注意，一定要用数学语言，精准的描述。

选择建模方法

其实，现在的数学领域已经将建模方法开发的比较成熟了，即使是我们平常碰到的NP问题，也都容易归结为NPC问题。

另一方面，基本的数学方法我们其实也都学的差不多了，即数学分析，微分方程，线性代数，概率论。

也就是说，其实就两个途径。

• 如果能自己想到好办法，可以自己先解决
• 想不到就去找文献，因为我们不是在做开拓性任务，如果真的做到开拓，那么就可以发论文了，还比个锤子的赛。

所以，由此可知，基础知识的储备+论文检索和利用的能力造就了最终的获奖。

而基础知识我们会在大一大二上学习，同时最好是自己看Mathematical modeling这种书来储备一些基本模型。

论文检索就得通过比赛去练。由此可知，不打比赛想上来就拿奖的，至少得需要有科研经历。否则就老实点找个普通比赛练手。

推导模型的数学表达式

• 用1步骤中的数学语言，把1中的符号写成2中的形式。

可能需要我们修改一些符号，我们可以通过赋值来解决，比如$f(t)=P$，其实就是换了个名字。

• 有可能需要补充假设，这个由找到的模型而定

新加进来的需要记好，个人感觉应该再补充到前面。

求解模型

• 模型已经套进去了，用3中已经数学化的方法去算结果即可。

注意推导不要有错误，似乎是一句废话。

• 运用计算机辅助代数，绘图，数值计算

这可以减少错误

回答问题

这一步是论文手和编程手的主打阶段，求解过程的可视化描述，以及论文的撰写。

这一步需要综合术语和通俗描述，达到雅俗共赏。

灵敏度分析

实际上，我们的计算有很多不确定度，假设并不一定稳定，比如价格的涨跌。

一般来说，我们都将这种不确定变量通过假设限定，也就是作为常量。但是因为他是可变的，所以实际上是超参数（hyperparameter）

灵敏度公式如下：

$$S(y,x)=\underset{△x\to 0}{\lim }\frac{\frac{△y}{y}}{\frac{△x}{x}}=\frac{dy}{dx}\cdot \frac{x}{y}$$

前面的式子是用来理解的：

直观来看，就是变化率的比值。比如x变1，y变3，那么就是三倍的变化率，相对来说影响比较剧烈。

他可以代表一个变量的变化对另一个变量的影响程度，用变化率去比较，更有代表性，可以消除绝对值差距过大带来的影响。

一般来说，我们的y主要是最优解对应的决策变量。x一般是超参数，就是不确定因子。

后面的式子是实际上计算的：

假设一个模型是可解的，假设只有x一个决策变量，r一个超参数，其他都看做定值。

那么一个目标函数y的最优解一定与x对应，可以通过最值条件直接求出x，x中带有r超参数，即x（r）。

也就是说，可以写出y和x的函数，也可以写出x和r，y和r的函数。

那么由此，无论是$\frac{dx}{dr}$还是$\frac{dy}{dr}$，$\frac{dy}{dx}$都是可以通过先得出一一对应的函数后求出来的。

也就是说，任意灵敏度都可以求出来。

关联变量灵敏度快速代换

很多时候，我们要求的两个变量，是有关系的，这个时候求灵敏度可以通过代换求解，比如P47页。

$$\begin{gathered} \because g=200c \\ \therefore S(x,g)=S(x,c) \end{gathered}$$

这个例子说明线性相关的两个变量灵敏度相同。你说为什么不是200倍？我们定义灵敏度为比例的比值就是为了消除这种绝对数值带来的影响，200倍大小不影响变化规律

又比如：
$$\because h=\frac{5}{c}$$
$$\therefore S(x,h)=\frac{dx}{dh}\cdot \frac{h}{x}$$
$$\begin{gathered} =\left(\frac{\frac{dx}{dc}}{\frac{dh}{dc}}\right)\left(\frac{\frac{5}{c}}{x}\right) \\ =-S(x,c) \end{gathered}$$

这个就是实现一个代数转化，这个在大物里面很常用。

实际上，从上面两个式子看出:

1. if $a\propto b$ then $S(x,a)=S(x,b)$
2. if $a\propto \frac{1}{b}$ then $S(x,a)=-S(x,b)$

灵敏度分析原则

不可能对每一个都分析，所以需要经验积累获得一种直觉。去判断出哪个超参数的不确定度最大。

灵敏性与稳健性

当超参数的灵敏性都很小，那么这个模型就是具有稳健性的，或者说鲁棒性。

如果灵敏度很大，按照我个人的理解，可能就要将他作为决策变量来去算了。至少不能放着不管，因为你的决断已经是不可靠的了。

多变量最优化

无约束最优化

所谓的无约束就是，仅仅对定义域做一个范围限制而无其他约束。

实际上就是找一个多维的最值即可。

具体执行

• 对n个决策变量求n个偏导，得到一个n元n行方程组
• 解方程即可
• 灵敏性分析，仍然满足前面的一一对应关系，即使是多了元。

理想情况应该是线性的方程组。要达成这个目标，f中就不应该有3次项，而不论二次项如何组合，求个导最后就只剩一次了，满足线性特征。

计算机辅助

Maple和Mathematical可以求解多元方程并且绘图，Matlab应该也行，但是还没有见过。

令我感觉意外的是，Maple上面的方程是直接接出来的，最后给出分数形式，并且可以添加超参数，有趣。

python的话，貌似得用梯度下降法，应该有现成的包。绘图用mpltoolkits的3D库。

拉格朗日乘子

原理

深入理解拉格朗日乘子法和KKT法

@知乎马同学 拉格朗日乘子法图形化理解

用我自己的理解来辅助一下记忆:

就是即使是没有约束条件，你也可以求出极值对吧，也就是求偏导后参数数量$\equiv$方程数量

理论上，只需要有一个函数包括所有参数，就可以求出最优。

那么我加了l个约束条件$h_k(x)=0$，怎么把这些统合到一个方程里呢？

$$F(x,\lambda )=f(x)-\sum _{k=1}^l{\lambda }_k{h}_k(x)$$

就是这么简单粗暴，引入一个约束，就引入一个参数(这个lambda就是拉格朗日乘子)，至于是多少，反正最后都能算出来，要记得前面强大的偏导方程组。

另一种写法：

是上面的梯度形式，都是向量，所以最后可以写成m个分式，形成和上面一样的方程组。

$$\nabla f={\lambda }_1\nabla g_1+\cdots +{\lambda }_n\nabla g_n$$

例：求$x^2+y^2+z^2=3,x=1$上$f(x,y,z)=x+2y+3z$的最大值点

先求f的偏导数，直接写成梯度算子：

$\nabla f=(1,2,3)$

然后写出两个约束条件的梯度算子：

$$\begin{gathered} \nabla g_1=(2x,2y,2z) \\ \nabla g_2=(1,0,0) \end{gathered}$$

g对应${\lambda }_1{\lambda }_2$的拉格朗日乘子，进行线性组合列方程：

$$\begin{gathered} 1=2x{\lambda }_1+{\lambda }_2 \\ 2=2y{\lambda }_1+0 \\ 3=2z{\lambda }_1+0 \end{gathered}$$

将这三个方程代入两个约束条件里，消去x，y，z会得到两个二元方程，求解即可。

不等式约束求解

将问题拆分为等号情况和不等号情况，比如$x\le 3$就分割为$x=3,x<3$两种情况。然后对每一种情况计算梯度。

我们一般是先求不等号情况：

不等号实际上就是规划出一个区域。实际上，两个三维的等号约束会变成一个二维的等号约束，比如两个面相交与一条曲线，但是区域是不会影响维度的。所以区域实际上没有强制的约束。

简言之，我们就当没有约束，使用无约束最优化。最后再判断一下极值点是否符合区域范围。

具体地，我们只需要假设没有区域，然后计算梯度为0的点，

如果这个点在区域内，那就保留下来，等待比较，或许这个极值会成为最值。

如果这个点在区域外，即区域内的$\nabla f\ne 0$，则说明在这个开区间内没有极值，自然不可能有最值，直接排除区域内情况。

然后就是等号约束情况：

有时候等号约束不是无限范围的，比如一个矩形的一条边，是有尽头的，所以仍然需要事后判断。

对线段使用拉格朗日乘子法，计算最优点，然后又会出现两种情况：

1. 极值点在线段内。保留下来即可。
2. 极值点在线段外，即$\nabla f\ne 0$，那么就只剩下端点可能是最值了，计算并保留下来等待比较

最后把我们所有的极值点比较找出最值点即可

同样，以上计算也可以在Maple或者Mathematical软件里实现。

灵敏性分析与影子价格

操作起来还是将a看做超参数，然后还是不变的一一对应。

也就是说S(x,a)求法不变。

所以我们后面主要谈S（y，a）

无约束情况下的灵敏度分析

S（y，a），这个和之前单变量的思路不太一样。

曾经的单变量，可以求出y和a的一一对应关系，但是当决策变量增加，一次性化简所有的决策变量不太可能，也不好操作，所以得出结论：

在多元情况下推荐使用链式法则解S（y，a）

$$\frac{df}{da}=\frac{\partial f}{\partial x_1}\cdot \frac{\partial x_1}{\partial a}+\frac{\partial f}{\partial x_2}\cdot \frac{\partial x_2}{\partial a}+\frac{\partial f}{\partial a}$$

因为无约束条件下求出来的最优解，可以确保
$$\frac{\partial y}{\partial x_1}\cdot \frac{\partial x_1}{\partial a}+\frac{\partial y}{\partial x_2}\cdot \frac{\partial x_2}{\partial a}\equiv 0$$

所以就有

$$\frac{dy}{da}=\frac{\partial y}{\partial a}$$

也就是，链式求导和形式上直接求“偏导”的结果相同。

拉格朗日情况下的灵敏度分析

详见P32底部

拉格朗日乘子的情况下，往往不能满足$\nabla f=0$的条件，所以对超参数求偏导，并不一定能保证
$\frac{\partial y}{\partial x_1}\cdot \frac{\partial x_1}{\partial a}+\frac{\partial y}{\partial x_2}\cdot \frac{\partial x_2}{\partial a}=0$

好消息是有一种特殊情况还是可以形成形式相同的结论的：

即，最优点落在边界的曲线约束上，且超参数不影响该曲线，则灵敏度形式同无约束优化

那么最优点的变化轨迹g 肯定是与$\nabla f$垂直的。

我们给出决策变量随着a变化的速度向量:

$$\frac{dx}{da}=(\frac{d{x}_1}{da},\frac{d{x}_2}{da})$$

这个向量肯定和$\nabla f=(\frac{\partial y}{\partial x_1},\frac{\partial y}{\partial x_2})$垂直，所以有：
$$\nabla f\cdot \frac{dx}{da}=(\frac{\partial y}{\partial x_1},\frac{\partial y}{\partial x_2})\cdot (\frac{d{x}_1}{da},\frac{d{x}_2}{da})=\frac{\partial y}{\partial x_1}\cdot \frac{\partial x_1}{\partial a}+\frac{\partial y}{\partial x_2}\cdot \frac{\partial x_2}{\partial a}=0$$

妙哉！

这种超参数一般对结果影响不大，这里给出一个直观理解

对于落在约束线上的最优解，如果超参数不影响约束线，那么无论超参数怎么变，最优解一定在约束线上。

约束线的特点就是，和$\nabla f$垂直，而这个方向是变化最快的。

那反过来，约束线是垂直方向，直观理解就是变化最慢的，所以这种情况求出来的S（y，a）应该很小，通常在10%-40%。

影子价格

约束条件的${\lambda }_i$，可以代表$g_i(x)=c_i$时的$S（f，c_i）$，证明在P37，直观理解，可以是在$c_i$变化的时候，其他约束不变，所以其他约束对f的变化没有影响，自然只有${\lambda }_i$有影响了，即$\frac{df}{d{c}_i}={\lambda }_i$

$\frac{dy}{dc}$的形式和我们前面的a不太一样，本质区别在于，c会影响曲线。

左半边就不会为0了，老实点求导算出来就好了。

我们注意到，左半边的值恰好等于拉格朗日乘子的值，这也不是偶然。因为我们的约束线现在在移动，而这两个移动的方向都是$\nabla$，所以梯度的比值就是变化率

$$\nabla f=\lambda \nabla g$$

$$\frac{dy}{dc}=\frac{\nabla f}{\nabla g}=\lambda$$

这个值叫影子价格，代表了生产规模扩大造成的利润扩大，如果开发新产品的利润不够影子价格，那还不如直接扩大原有规模。

问题又来了，右半边为什么就0了？

那是因为，c是在约束条件里的，并且没有和x同时出现在y中，所以$\frac{\partial y}{\partial c}=0$

深入理解链式法则

参考

链式法则的理解可以到吴恩达深度学习的一节课里看，讲的挺明白的，我这里从形式上辅助记忆

从这个视频开始看完计算图

图形化理解链式法则

基础写法

基础的链式法则比较简单，对于多个中间变量和多个基础变量的原则如下

1. 有n基础变量，就要求n个偏导，形成n个方程
2. 对其中一个方程，假设其对应的基础变量为为a，然后假设有$x_1,x_2$两个中间变量。那么链式法则就可以写成下面的形式
   $$\frac{\partial f}{\partial a}=\frac{\partial f}{\partial x_1}\cdot \frac{\partial x_1}{\partial a}+\frac{\partial f}{\partial x_2}\cdot \frac{\partial x_2}{\partial a}$$

特殊情况

一般来说，公式里面不会出现基础变量，但是有的公式偏偏就把基础变量和中间变量写在了一起，即$f(x_1,x_2,a)$，这种求导需要多加一项。

即
$$\frac{\partial f}{\partial a}=\frac{\partial f}{\partial x_1}\cdot \frac{\partial x_1}{\partial a}+\frac{\partial f}{\partial x_2}\cdot \frac{\partial x_2}{\partial a}+\frac{\partial f}{\partial a}\cdot \frac{\partial a}{\partial a}$$

因为$\frac{\partial a}{\partial a}\equiv 1$，所以式子就是

$$\frac{\partial f}{\partial a}=\frac{\partial f}{\partial x_1}\cdot \frac{\partial x_1}{\partial a}+\frac{\partial f}{\partial x_2}\cdot \frac{\partial x_2}{\partial a}+\frac{\partial f}{\partial a}$$

这么做合理吗？合理，即使你用这个办法去计算没有基础变量的情况，也满足。

所以说后一种应该是前一种的泛化。

最优化计算方法

前面的两个，都是可以求解出来的，那两个也都是线性的方程组和约束，而现实往往是非线性的，要么就是计算非常繁杂，要么就是干脆不能代数求解，这时就应该通过数值近似求解，也就是数值分析，嘿嘿，课是好课，就是学校考核有问题罢了。

单变量最优化

流程

先用粗略的全局图像法：

图像法，先用大步长画图，选定区间，逐步迭代精细化放大，一般两三次就足够说明问题了。

然后分区间使用局部快速迭代解法：

然后先求导出导数公式再用牛顿法求解导数零点，或者直接梯度下降。

灵敏度分析

因为我们使用了数值方法求解，所以灵敏度不能用导数直接求解了，我们回归最初始的定义，就是一个差商公式：

$$\underset{△\to 0}{\lim }\frac{\frac{△x}{x}}{\frac{△r}{r}}$$

我们直接给r 来一个$ϵ$ 扰动，然后用差商计算即可。

其他计算可以使用这个灵敏度代换计算，不一定要直接算。

多变量最优化

变量增加后，图像法不再适用，求导也更加麻烦。有约束情况下，几何形状更加复杂，求解也更加困难。

先进行随机搜索防止局部最优，然后进行多变量牛顿法精确求解。

随机搜索算法

很朴素，很粗略的全局搜索。

有一种方法是进行网格搜索，选定步长进行一个伪枚举，然后逐步缩小范围，但是空间一大，这个就不行了。

随机搜索就是每次在解空间里随机生成一个点，然后判断大小来更新最值，迭代N次得出最终结果。

N越大，越准确，精确程度据说是和均匀分点差不多（理论上随机$\approx$均匀），极端情况下，如果随机的点几乎将所有解空间覆盖，那一定是最优解（然鹅这就是暴力遍历了）

与随机搜索匹配的，还有粗化解空间，缩小解空间的技巧。

比如解空间为[1,2,3,4,5,6,7]，那么我们可以粗化为[2,4,6]，用这个来代替周围的领域，简化搜索。

大致确定解空间后也可以缩小范围继续搜索，不过没有太大意义，我们该采取更精确的方法。

多变量牛顿法

只是将单变量扩展到了多变量，进行向量化，形式不变。

具体拆解一下向量化：

$x_{n\times 1}=\left[\begin{array}{c}{x}_1\\ x_2\\ \cdots \\ x_n\end{array}\right]$,$F(x{)}_{n\times 1}=\left[\begin{array}{c}{f}_1(x)\\ f_2(x)\\ \cdots \\ f_n(x)\end{array}\right]=\left[\begin{array}{c}{f}_1(x_1,x_2,\cdots ,x_n)\\ f_2(x_1,x_2,\cdots ,x_n)\\ \cdots \\ f_n(x_1,x_2,\cdots ,x_n)\end{array}\right]$

$A_{m\times n}=\left[\begin{array}{c}\nabla f_1\\ \nabla f_2\\ \cdots \\ \nabla f_n\end{array}\right]=\left[\begin{array}{ccc}\frac{\partial f_1}{\partial x_1}& \cdots & \frac{\partial f_1}{\partial x_n}\\ ⋮& & ⋮\\ \frac{\partial f_n}{\partial x_1}& \cdots & \frac{\partial f_n}{\partial x_n}\end{array}\right]$

线性近似公式为：

$$F_{n\times 1}(x)\approx F_{n\times 1}(x(0))+A_{n\times n}(x-x(0){)}_{n\times 1}$$

令F=0，得出计算公式为：

$$\begin{gathered} x=x-A^{-1}F(x) \\ {x}_{n\times 1}=x_{n\times 1}-A_{n\times n}^{-1}{F}_{n\times 1}(x(0)) \end{gathered}$$

当A是1*1大小的时候，取逆就是倒数，所以是向下兼容单变量的。必然是方阵，因为n个未知数有n个偏导方程。

至于求逆就交给计算机了。

线性规划

有约束的多变量最优化往往是最难求解的，特殊问题有特殊求法，但是目前还缺乏一般性的方法。

这种问题的一个简化就是线性问题，即：

1. 目标函数是线性
2. 约束函数是线性

求解流程

单纯形法理解

合理性：

1. 一般来说，我们线性规划都是给不等式的，但是实际上可行域内的$\nabla f$不会为0（存疑），所以一定是在边界，同样是因为梯度不为零，所以边界会出现在顶点处。
2. 单纯形法实际上就是逐个试点，求最小的顶点。

线性规划标准模型为：

1. $x=(x_1,\cdots ,x_n)$
2. $g_i(x)\le c_i$为约束，假设有m个约束
3. $x_i\ge 0$，这是n个固有约束。
4. 则总共上有m+n个约束。

因为将不等式转化一下，变成单纯形法的标准模型：

1. 将不等号变为等号
2. 在左边补充一个松弛变量，m个约束各一个总共补充m个，当松弛变量为0时，i约束为关键约束。

大致求解过程：

1. 从m+n个变量中，选出n个变量置0，剩下m个变量m个方程，可以线性求解，这样就诞生了一个顶点。
2. 用$x_i\ge 0$来进行判断，全部非负则为可行解。
3. 然后逐个顶点测试

优化：

1. 从上面可以看出，复杂度为$C_{m+n}^n$，复杂度是很高的，所以实际上单纯形法只是计算一个顶点子集。

我们的操作：

1. 实际中，我们不需要去计算，我们只需要把线性规划写成标准的模式即可，然后由计算机实现具体求解。

灵敏性估计

操作比较简单，就是在程序里改一下数据，做一个差分计算即可。

走到这个复杂度，我么已经不能够代数求灵敏性了，顶多是进行一个转换。

以P61对燕麦的灵敏性分析进行理解：

1. 这里加深了对超参数的理解。本质上来说，在约束条件里的超参数，改变的是解空间，而在目标函数里的超参数，改变的是梯度分布。解空间和梯度分布的综合决定了最优解的位置，综合决定了影子价格！
2. 基于1，会产生一个特殊的情况，即梯度和两顶点连线垂直，因为线性规划的边界是线性的（可以理解为直线），所以在两点上移动不会影响结果，即条线上的任一点其实都是最优解，只不过程序默认选择端点罢了。
3. 这种情况，结果就要对两个端点进行t和1-t的线性组合，再结合其他实际情况来进行调整，来同时达到另一个目标函数的效益最大化。

电子表格分析后的体会：

1. 书里说使用Quattro Pro来分析灵敏度好，我感觉就那样，如果要用的话也建议Excel了
2. 影子价格体现了在约束下解的性质，如果影子价格突变，那么解的性质就不一样了，直接体现就是决策变量变了，代表了策略的改变，空间上说就是最优解换顶点了。反过来说，还是有一定缓冲范围的，这个缓冲范围就是增加量/减少量

稳定性估计：

1. 对决策结果想进行进一步的分析，比如想知道为什么不进行彻底利用，如果不留下来会怎样，这种问题就可以自己改一下约束来看看结果
2. 如果结果差距不大，可以制定出多个决策，来适应不同情况，毕竟，真实的约束远比模型多

离散最优化

有一些问题是离散的，可行域比较复杂，经常用图或者树结构来描述，但是目前也仍然没有普遍有效的方法。

我个人感觉，最难受的地方在于没法求导了，或者是只能粗略的进行差分。

线性规划（LP问题）失去连续性后细化为整数规划（IP问题），整数规划也属于基本问题，同时和线性规划比较相似，便于对比学习。

整数规划-基本流程

分支定界法解释

分支定界核心思想：

1. 解空间实际上可以被不断二分直到每一种离散的可能，形成一个二叉树搜索空间。但是这个空间实在是太大了，需要剪枝。
2. 剪枝的原则就是，解除当前解空间的离散限制，假定为连续的，转化为LP松弛问题，用线性规划的单纯形法去求解，如果连续情况的最大值都低于我们当前的最大值，那已经没希望了，剪枝，不再进行递归二分。
3. 另一种提前终止的情况就是走到了最大深度，已经满足所有的离散约束，这就是一个可行解，我们进行比较，如果更优，就更新。
4. 最后就可以找出最优解，通过剪枝，我们省去了很多没有希望的搜索。但是即使是剪枝，我们也进行了很多次线性规划，所以IP问题耗时挺严重。

二值整数规划模型（BIP）：

当离散值为0，1时，已有很多快速的算法。

灵敏性分析

灵敏性分析很麻烦，原因如下：

1. 我们要重新计算，计算IP问题本身就耗时
2. 最优解不一定（大概率）不出现在边缘上，结果就是不光滑，没办法求梯度之类的，没有影子价格，不过可以用边际效益代替。
3. 离散的结果就是，要么不变，要么突变，0和1就是有和无的区别，1和2就是一倍的差别。

对可行域离散化的理解：

1. 实际上，如果切分得非常细，就是精确到0.01的级别，你会发现离散情况的最优解是和连续的最优解几乎一样的
2. 所以本质上说，连续就是粒度无限细的离散情况。
3. 当粒度变粗，逐渐离散，离散最优解和连续最优解有大偏差的几率就会显著增加。