《数学建模算法与应用》司守奎著 阅读笔记。
一、线性规划
线性规划:在一组线性约束条件的限制下,求一线性目标函数的最大或最小值
一般标准型:
Matlab标准型:
Matlab求解线性规划的命令为:
[x,fval] = linprog(f,A,b,Aeq,beq,lb,ub)
可转化为线性规划的问题:
取,
![]()
,
![]()
:
二、整数规划
整数规划分为两大类:1)变量全限制为整数,纯整数规划;2)变量部分限制为整数,混合整数规划。
注:整数规划的最优解不能按照实数最优解简单取整而获得。
2.1 0-1型整数规划
0-1型整数规划是将某些变量
![]()
限制为二进制变量的规划问题。在实际问题中,引入0-1变量可以把需要
分类讨论的数学规划问题
统一为一个规划问题。下面举例:
1)相互排斥的约束条件规划问题:
如果有
![]()
个互相排斥(需要分类讨论)的约束条件
由于
![]()
个约束条件互相排斥,故只有一个会起作用,引入
![]()
个0-1变量来避免分类讨论:
引入一个充分大的常数
![]()
,可将
![]()
个互相排斥(需要分类讨论)的约束条件,统一为
![]()
个约束条件:
由(7)式知,当
![]()
时,由于M充分大,则不等式恒成立,约束条件失效;由(8)式可知,仅有一个
![]()
可以取
![]()
.
2)固定费用问题和指派问题
见书P13,P14 ,不太好统一形式,不过多赘述了。。
2.2 蒙特卡洛法(随机取样法)
蒙特卡洛法是基于大量事件的统计结果来实现一些确定性问题的计算方法。
下面举一个实际编程案例:
例:
![]()
与
![]()
轴在第一象限围成的一个曲边三角形。用随机实验法求得该三角形的面积。
解:随机实验的思想为,在
![]()
的矩阵区域内,按均匀分布取样若干大量的随机点,统计落入计算区域的随机点个数,求得频数,曲边三角形的面积则为矩阵面积乘频数。
clc,clear
x = unifrnd(0,12,[1,10000000])
y = unifrnd(0,9,[1,10000000])
pinshu = sum(y<x.^2 & x<=3) + sum(y<12-x & x>=3);
area_appr = 12 * 9 * pinshu/10^7
非线性整数规划问题:
非线性规划问题尚未具备一种通用成熟且准确的求解方法,但由于加入了整数限制,则非线性规划问题的解的个数也变为了有限个,有了蒙特卡洛法后,我们则可以用枚举法求出满意解。
具体编程案例见书P15。
2.3 整数线性规划的计算机求解
Matlab下混合整数线性规划的标准型为:
其中
![]()
为整数变量的索引。
Matlab求解整数线性规划的命令:
[x,fval] = intlinprog(f,intcon,A,b,Aeq,beq,lb,ub)
三、非线性规划
非线性规划:目标函数或约束条件包含非线性函数的规划问题。
非线性规划的数学模型一般形式:
其中:
![]()
为模型的决策变量;
![]()
为目标函数;
![]()
为等式约束;
![]()
为不等式约束;且
![]()
中包含非线性函数。
注:若线性规划问题的最优解存在,则此最优解只能在可行域的边界上达到;而非线性规划的最优解可能在可行域的任意一点达到。
Matlab下非线性规划的标准形式:
其中
![]()
为目标函数;
![]()
为非线性向量函数。
Matlab求解此标准型非线性规划问题的命令为:
[x,fval] = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)
其中,fun为用M文件定义的目标函数
![]()
; nonlcon是用M文件定义的非线性向量函数
![]()
; options 定义了优化参数,可用Matlab默认设置;其余变量若无则用"[]"代替
3.1 无约束极值问题
无约束极值问题:即没有约束条件,仅求解目标函数极值的优化问题。
其数学模型可简单的表示为:
在Matlab中有两个专用函数求解该极小值问题,fminunc 和 fminsearch。
Matlab中 fminunc的基本命令:
[x,fval] = fminunc(fun,x0,options)
其中:fun是用M文件定义的函数。当fun只有一个返回值时,是函数值
![]()
; 当 fun 有两个返回值时,它的第二个返回值是
![]()
的梯度向量;当fun 有三个返回值时,它的第三个返回值是
![]()
的二阶导数矩阵(Hessian阵)。x0是
![]()
的初始值。
Matlab中 fminsearch的基本命令:
[x,fval] = fminsearch(fun,x0,options)
注:fminsearch只能求给定初始值附件的一个极小值点。
3.2 约束极值问题
约束极值问题:带约束条件的极值问题,又称规划问题。
3.2.1 二次规划
二次规划问题:目标函数为自变量
![]()
的二次函数,约束条件是线性的规划问题。
Matlab中二次规划的数学模型标准形式:
其中:
![]()
为实对称矩阵。
Matlab中求解二次规划的命令是:
[x,fval] = quadprog(H,f,A,b,Aeq,beq,lb,ub,x0,options)
(具体细节可参考Matlab命令窗口中运行 help quadprog 后的“帮助”文档)
3.2.2 罚函数法
罚函数法: 将带约束的非线性规划问题,转化为一系列无约束极值问题。
(类似于拉格朗日乘数法,但也区别于拉格朗日乘数法)
罚函数法求解非线性规划问题的思想是,利用问题中的约束函数作出适当的罚函数,由此构造出带参数的增广目标函数,把问题转化为无约束非线性规划问题。
外罚函数法的例子:
取一个充分大的数
![]()
,构造函数
不难看出,为了让增广目标函数
![]()
取得极值,则
![]()
,
![]()
三个条件必须满足。
3.3 Matlab求解约束极值问题
1.fminbnd 函数
求单变量非线性函数在区间上的极小值:
Matlab的命令为:
[x,fval] = fminbnd(fun,x1,x2,options)
该命令中,fun为用M文件定义的单变量数学函数。
2. fseminf 函数
求含附加变量标量函数约束的非线性规划问题:
其形式与非线性规划标准形式(11)式类似,而其中
![]()
为标量函数,
![]()
为附加变量。
上述问题的Matlab求解命令为:
[x,fval] = fseminf(fun,x0,ntheta,seminfcon,A,b,Aeq,beq,lb,ub)
其中,ntheta是半无穷约束
![]()
的个数;函数seminfcon用于定义非线性不等式约束
![]()
、非线性等式约束
![]()
和半无穷约束
![]()
的函数,函数seminfcon有两个输入参量x 和s , s 是推荐的取样步长,可不使用。
3. fminimax 函数
求解
Matlab的求解命令为:
[x,fval] = fminimax(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)
3.4 Matlab优化工具箱的图形界面解法
Matlab优化工具箱的optimtool命令提供了优化问题的用户图像界面解法。
具体操作步骤见书P32
四、 图与网络模型及其方法
图论中所谓的"图"是指某类具体事物的这些事物之间的联系。如果用点表示这些具体事物,用连接两点的线段(直的或曲的)表示两个事物的特定联系,就得到了描述这个"图"的几何形象。
4.1 图的基本概念与数据结构
4.1.1 基本概念
平面上有
![]()
个点,用直线或曲线将这些点连接起来,不考虑点的位置与连线曲直长短,这样形成的一个关系结构就是一个图,记作
![]()
,其中
![]()
是上述点为元素的顶点集,
![]()
是以上述连线为元素的边集。
有向图: 各边都包含方向的图结构。
无向图:各边都不含方向的图结构。
混合图:有的边有方向,有的边无方向的图结构。
简单图:任两顶点间最多有一条边,且每条边的两个端点皆不重合的图结构。
完全图:如果图的两顶点间有边相连则称两顶点相邻,每一对顶点都相邻的图称为完全图,否则为非完全图。
二分图:若顶点集
![]()
(这里
![]()
表示顶点集
![]()
中元素的个数,即两个顶点集的元素个数均不为零),且
![]()
中无相邻顶点对,
![]()
中亦然,则称图
![]()
是二分图;特别地,若任意
![]()
都与
![]()
中每个点相邻,则称图
![]()
是完全二分图,记作
![]()
。(
就是一个图的顶点能分为两组,组内点不能相邻,组间点能相邻)
度: 设点
![]()
是边
![]()
的端点,则称
![]()
与
![]()
相关联,与
顶点
相关联的边数称为该顶点的度,记为
![]()
。度为奇数的点称为奇顶点,度为偶数的顶点称为偶顶点,不难想到
![]()
,即
所有顶点的度数之和是边数的2倍,由此可知
奇顶点的总数是偶数。
道路:设
![]()
,称
![]()
是图
![]()
的一条道路,
![]()
为路长,
![]()
为起点,
![]()
为终点。
迹、轨道、圈:各边相异的道路称为迹;各顶点相异的道路称为轨道,可记作
![]()
;
起点与终点重合的
道路称为
回路;
起点与终点重合的
轨道称为
圈,即
![]()
中的
![]()
。
连通图、Hamilton轨、Hamilton圈、Hamilton图:图中任两顶点之间都存在道路的图,称为连通图。图中含有所有顶点的轨道称为Hamilton轨;闭合的Hamilton轨称为Hamilton圈;含有Hamilton圈的图称为Hamilton图。
距离:称两顶点
![]()
分别为
起点和终点的最短轨道之长为顶点
![]()
的
距离;在
完全二分图
中,
中的两顶点之间的距离为偶数,
与
中的顶点之间的距离为奇数。
赋权图:每条边上都有一个(或多个)实数对应的图,这个(些)实数称为这条边的权。
4.1.2 图与网络的数据结构(表示法)
为了在计算机上实现网络优化的算法,首先必须有一种方法(即数据结构)在计算机上描述图与网络。
计算机上图与网络的表示方法主要有两种:邻接矩阵表示法和稀疏矩阵表示法。
设
![]()
是一个简单无向图,顶点集合
![]()
,边集
![]()
,记
![]()
,
![]()
。
1.邻接矩阵表示法
邻接矩阵是表示顶点之间相邻关系的矩阵,邻接矩阵记为
![]()
。
![]()
为赋权图时:
![]()
为非赋权图时:
缺点: 当图的边数
远小于顶点数
时,邻接矩阵表示法会造成很大的空间浪费。2.稀疏矩阵表示法
稀疏矩阵是指矩阵中零元素很多,非零元素很少的矩阵。
稀疏矩阵只需存储非零元素的位置,即(非零元素的行地址,非零元素的列地址,非零元素的值)
在有向图中,Matlab可直接使用sparse命令将邻接矩阵转化为稀疏矩阵表示方式。
在无向图中,Matlab只需存储邻接矩阵的下三角元素,因为邻接矩阵为对称矩阵。
稀疏矩阵只是一种存储形式,使用时通过full命令变回普通矩阵。
4.2 最短路径问题
4.2.1 指定顶点到其他各顶点之间的最短路径
问题描述:给出一个赋权图网络结构,每边上的权代表两顶点间的距离,求出指定起点到其他各顶点的最短路径。
迪克斯特拉(Dijkstra)算法
原书的算法描述不太好懂,用更通俗的描述记录一下。。。
Dijkstra算法是一种贪心的策略,简单来说就是从起点开始,找到与其路径最短的点并存储进已找出最短路径的顶点集合
![]()
,再用此顶点集合的点寻找新的路径(路径只能经过顶点集合
![]()
中的点)并更新路径长度,寻找下一个路径最短的点,以此类推直到所有顶点都找出最短路径。
具体过程先看如下示例: (来源:Dijkstra算法原理)
来源:Dijkstra算法原理求
![]()
中以顶点
![]()
为起点到其他各点的最短路径距离:
来源:Dijkstra算法原理现在再看原书中的算法描述:
设图
![]()
中,起点为
![]()
,
![]()
表示点
![]()
到顶点
![]()
的
当前最短路径,第
![]()
次迭代后已找到最短路径的顶点集合为
![]()
。
初始时,
(1) 对每个
![]()
(指集合
![]()
中除去
![]()
的元素),用下式
更新所有
![]()
的值,并将其中值最小
![]()
,即最短路径的顶点
![]()
记为
![]()
存入顶点集合
![]()
,得到新的顶点集合
![]()
。
(2) 令
![]()
, 循环操作步骤(1)
(3) 直到
![]()
时,遍历结束,起点
![]()
到各顶点
![]()
的最短路径为存入
![]()
之前最后一次计算求得的
![]()
值。
为避免重复计算,算法过程中采用标号法,当
![]()
进入顶点集
![]()
之前标签为
![]()
标号,当
![]()
进入顶点集
![]()
之后标签为
![]()
标号。当每个顶点都标号为
![]()
时,算法结束。
该算法在Matlab中的具体实现为:
从起点sb到终点db的Dijkstra算法程序:
function [mydistance,mypath] = mydijkstra(a,sb,db);
% 输入:a——邻接矩阵;a(i,j)——i到j之间的距离,可以是有向的
% sb——起点的标号,db——终点的标号
% 输出:mydistance——最短路径的距离,mypath——最短路径
n = size(a,1);visited(1:n) = 0; % visited==0指标号T,visited==1指标号P
distance(1:n) = inf; distance(sb) = 0; %起点到各顶点距离的初始化
visited(sb) = 1; u = sb; % u为最新的P标号顶点
parent(1:n) = 0 % 前驱顶点的初始化,即顶点集S
for i = 1:n-1
id = find(visited == 0); % 标号T的顶点
for v = id
if a(u,v) + distance(u) < distance(v) % 即min{ l(v), l(u)+w(uv)
distance(v) = distance(u) + a(u,v) % 更新最短距离
parent(v) = u ;
end
end
temp = distance;
temp(visited == 1) = inf; % 标号P的点临时距离变量设置为无穷大(防止其重复计算)
[t,u] = min(temp); % 求出最短距离的顶点,记为新的u
visited(u) = 1; % 将新的u标记 为P点
end
mypath = [];
if parent(db) ~= 0 % 若存在路径
t = db; mypath = [db];
while t ~= sb
P = parent(t);
mypath = [P mypath];
t = P;
end
end
mydistance = distance(db);
4.2.2 每对顶点之间的最短路径
Dijkstra算法可计算指定顶点到其他各顶点的最短路径,若要计算赋权图中所有顶点两两成对的最短路径,则需要对每个顶点循环调用Dijkstra算法求出以该点为起点到其余各点的最短路径。但这样的方法需要执行n-1次,时间复杂度为
![]()
,稍微有点笨拙。
另一种解决每对顶点之间最短路径的算法为——Floyd算法。
对于赋权图
![]()
,初始时定义邻接矩阵:
式中,
对于无向图,
![]()
是对称矩阵,
![]()
。
Floyd算法的基本思想是递推产生一个矩阵序列
![]()
,其中矩阵
![]()
的第
![]()
行第
![]()
列元素
![]()
表示从顶点
![]()
到顶点
![]()
的路径上所经过的顶点序号不大于
![]()
的最短路径长度。
计算时用迭代公式
![]()
是迭代次数,
最后,当
![]()
时,
![]()
即是各顶点之间的最短通路值。
Matlab通用解法,起点sb到终点db通用的Floyd算法程序如下:
function[dist, mypath] = myfloyd(a, sb, db)
% 输入:a——邻接矩阵; 元素a(i,j)——顶点 i 到 j 之间的直达距离,可以是有向的
% sb——起点的标号; db——终点的标号
% 输出:dist——最短路的距离; % mypath——最短路的路径
n = size(a,1); path = zeros(n);
for k=1:n
for i=1:n
for j=1:n
if a(i,j) > a(i,k) + a(k,j)
a(i,j) = a(i,k) + a(k,j);
path(i,j) = k;
end
end
end
end
dist = a(sb,db);
parent = path(sb,:); % 从起点 sb 到终点 db 的最短路上各顶点的前驱顶点
parent(parent == 0) = sb; % path 中的分量为0,表示该顶点的前驱是起点
mypath = db; t = db;
while t ~= sb
p = parent(t);mypath = [p,mypath];
t = p;
end
4.3 最小生成树问题
树:连通的无圈图称为树,记作
![]()
,图
![]()
是树的充要条件为
![]()
中任意两顶点之间有且仅有一条轨道。
最小生成树:在连通赋权图上求权值和最小的生成树,称为最小生成树。
最小生成树的两种常用求解算法:prim算法、Kruskal算法。
4.3.1 prim算法
prim算法的基本思路为贪心的寻找当前最小权值的边,递归地将所有顶点包含。
prim算法定义两个集合P和Q,其中P存放G中最小生成树的顶点,Q存放G中最小生成树的边。初始时,
![]()
(
![]()
为起点),
![]()
(空集)。从
![]()
和
![]()
中分别选择起点和顶点找出最小权值的边记作
![]()
,将该顶点
![]()
加入集合P, 边
![]()
加入集合Q。重复循环直到
![]()
时,最小生成树构造完毕。集合
![]()
包含最小生成树的所有边。
prim算法流程:
①
② while
![]()
:
![]()
;
![]()
.
end
4.3.2 Kruskal算法
Kruskal算法流程如下:
① 选择
![]()
,且
![]()
是权值最小的边。
② 若
![]()
已选好,从
![]()
中选择
![]()
,且
![]()
需满足:
③ 直到选得
![]()
为止。
4.4 最大流问题
网络: 给一个有向图
![]()
,其中
![]()
为弧集,弧包含一个发点
![]()
,一个收点
![]()
,其余点为中间点,弧
![]()
, 每条弧对应一个容量值
![]()
(简称
![]()
),称为弧的容量。通常把这样的有向图
![]()
叫做网络,记为
![]()
,其中
![]()
.
流: 指定义在弧集合
![]()
上的一个函数
![]()
, 并称
![]()
为弧
![]()
上的流量。
可行流:可行流需要满足两种约束条件:容量限制条件和流量平衡条件
① 容量限制条件: 对每一弧
![]()
。
②流量平衡条件:
中间点
![]()
(既有流入也有流出):流出量=流入量,即对于每个
![]()
,有
发点
![]()
:流出量=可行流的流量
![]()
,有:
收点
![]()
: 流入量= - 可行流的流量
![]()
,有:
则,最大流问题可以写为如下的线性规划模型:
4.4.1 最小费用最大流问题
给定网络
![]()
,每一弧
![]()
上,除了已给容量
![]()
外,还给了一个单位流量的费用
![]()
(简记为
![]()
)。而所谓最小费用最大流问题,则是希望求得一个从发点
![]()
到收点
![]()
的最大流,且使最大流的总输送费用
![]()
取最小值。因此,最小费用最大流问题,相当于两个线性规划模型,即在求解式(26)后,再求出相应最大流所对应的最小费用,即求解以下线性规划问题。
其中:
![]()
为线性规划模型(26)求解所得的最大流的流量。
4.4.2 Matlab工具箱求解最大流问题
最大流问题和最小费用最大流问题均可以通过Matlab工具箱进行求解。
Matlab图论工具箱的命令见下表。
命令名 | 功能 |
---|
graphallshortestpaths | 求图中所有顶点对之间的最短距离 |
graphconncomp | 找无向图的连通分支,或有向图的强(弱)连通分支 |
graphisdag | 测试有向图是否含有圈,不含圈返回1,否则返回0 |
graphisomorphism | 确定两个图是否同构,同构返回1,否则返回0 |
graphisspantree | 确定一个图是否是生成树,是返回1,否则返回0 |
graphmaxflow | 计算有向图的最大流 |
graphminspantree | 在图中找最小生成树 |
graphpred2path | 把前驱顶点序列变成路径的顶点序列 |
graphshortestpath | 求图中指定一对顶点间的最短距离和最短路径 |
graphtopoorder | 执行有向无圈图的拓扑排序 |
graphtraverse | 求从一顶点出发,所能遍历图中的顶点 |