数学建模算法大全


数 学 建 模 主 编 司守奎 副主编 徐珂文 李日华 海军航空工程学院 ﹒i﹒ 目 录 第一章 线性规划……………………………………………………………………………1 § 1 线性规划……………………………………………………………………………1 1.1 线性规划的实例与定义 ……………………………………………………1 1.2 线性规划的 Matlab 标准形式…………………………………………………1 1.3 线性规划问题的解的概念…………………………………………………… 2 1.4 线性规划的图解法…………………………………………………………… 2 1.5 求解线性规划的 Matlab 解法…………………………………………………3 1.6 可以转化为线性规划的问题………………………………………………… 4 § 2 运输问题………………………………………………………………………… 4 § 3 指派问题………………………………………………………………………… 5 § 4 对偶理论与灵敏度分析………………………………………………………… 7 习题一……………………………………………………………………………… 9 第二章 整数规划………………………………………………………………………… 12 § 1 概论……………………………………………………………………………12 § 2 分枝定界法………………………………………………………………………12 § 3 10 − 整数规划………………………………………………………………… 14 3.1 引入 10 − 变量的实际问题………………………………………………… 15 3.2 10 − 整数规划解法之一………………………………………………… 16 § 4 蒙特卡洛法(随即取样法)……………………………………………………17 §5 整数规划的计算机解法…………………………………………………………18 习题二………………………………………………………………………………19 第三章 非线性规划………………………………………………………………………20 § 1 非线性规划………………………………………………………………………20 1.1 非线性规划实例与定义 ……………………………………………………20 1.2 线性规划与非线性规划的区别 ……………………………………………21 1.3 非线性规划的 Matlab 解法…………………………………………………21 1.4 求解非线性规划的基本迭代格式 …………………………………………22 1.5 凸函数、凸规划…………………………………………………………… 23 § 2 无约束问题……………………………………………………………………… 23 2.1 一维搜索方法……………………………………………………………… 23 2.2 二次插值法………………………………………………………………… 26 2.3 无约束极值问题的解法 ……………………………………………………26 2.4 Matlab 求无约束极值问题……………………………………………… 32 § 3 约束极值问题…………………………………………………………………32 3.1 二次规划…………………………………………………………………… 33 3.2 罚函数法…………………………………………………………………… 34 ﹒ii﹒ 3.3 Matlab求约束极值问题…………………………………………………………35 § 4 飞行管理问题…………………………………………………………………38 习题三………………………………………………………………………………………39 第四章 动态规划…………………………………………………………………………40 § 1 引言………………………………………………………………………………40 § 2 基本概念,基本方程和计算方法………………………………………………41 § 3 逆序解法的计算框图……………………………………………………………43 § 4 动态规划与静态规划的关系……………………………………………………44 § 5 若干典型问题的动态规划模型…………………………………………………46 § 6 具体的应用实例…………………………………………………………………47 习题四…………………………………………………………………………………50 第五章 图与网络模型及方法…………………………………………………………… 51 § 1 概论………………………………………………………………………………51 § 2 图与网络的基本概念……………………………………………………………52 § 3 应用—最短路问题………………………………………………………………58 § 4 树…………………………………………………………………………………60 § 5 匹配问题…………………………………………………………………………63 § 6 Euler图和 Hamilton 图………………………………………………………… 64 § 7 最大流问题………………………………………………………………………68 § 8 最小费用流及其求法……………………………………………………………73 习题五……………………………………………………………………………… 74 第六章 排队论模型……………………………………………………………………… 76 § 1 基本概念…………………………………………………………………………76 § 2 输入过程与服务时间的分布……………………………………………………78 § 3 标准的 M/M/1 模型………………………………………………………………81 § 4 产生给定分布的随机数的方法…………………………………………………82 § 5 排队模型的计算机模拟…………………………………………………………83 习题六………………………………………………………………………………… 86 第七章 对策论…………………………………………………………………………… 87 § 1 引言………………………………………………………………………………87 § 2 对策问题…………………………………………………………………………87 § 3 零和对策的混合策略……………………………………………………………90 § 4 零和对策的线性规划解法………………………………………………………92 习题七……………………………………………………………………………… 95 第八章 层次分析法……………………………………………………………………… 96 § 1 层次分析法的基本原理与步骤……………………………………………… 96 § 2 层次分析法的应用……………………………………………………………100 习题八……………………………………………………………………………… 103 第九章 插值与拟合……………………………………………………………………… 104 ﹒iii﹒ § 1 插值方法………………………………………………………………………104 1.1 拉格朗日多项式插值………………………………………………………104 1.2 牛顿插值………………………………………………………………………106 1.3 分段线性插值…………………………………………………………………108 1.4 埃尔米特(Hermite)插值………………………………………………………109 1.5 样条插值………………………………………………………………………110 1.6 二维插值………………………………………………………………………113 § 2 曲线拟合的线性最小二乘法……………………………………………………114 2.1 线性最小二乘法………………………………………………………………114 2.2 最小二乘法的 Matlab 实现…………………………………………………… 116 § 3 最小二乘优化……………………………………………………………………117 § 4 曲线拟合与函数逼近…………………………………………………………119 习题九…………………………………………………………………………………120 第十章 数据的统计描述和分析……………………………………………………122 § 1 统计的基本概念………………………………………………………………122 § 2 参数估计………………………………………………………………………128 § 3 假设检验……………………………………………………………………129 习题十…………………………………………………………………………………133 第十一章 方差分析…………………………………………………………………134 § 1 单因素方差分析………………………………………………………………134 § 2 双因素方差分析………………………………………………………………138 习题十一………………………………………………………………………………140 第十二章 回归分析…………………………………………………………………141 § 1 多元线性回归…………………………………………………………………141 § 2 非线性回归和逐步回归………………………………………………………148 习题十二………………………………………………………………………………151 第十三章 微分方程建模……………………………………………………………153 § 1 发射卫星为什么用三级火箭…………………………………………………153 § 2 人口模型………………………………………………………………………158 § 3 战争模型………………………………………………………………………160 习题十三………………………………………………………………………………165 第十四章 稳定状态模型……………………………………………………………167 § 1 微分方程稳定性理论简介……………………………………………………167 § 2 再生资源的管理和开发………………………………………………………169 § 3 Volterra 模型………………………………………………………………… 174 习题十四………………………………………………………………………………178 第十五章 常微分方程的解法………………………………………………………179 § 1 常微分方程的离散化…………………………………………………………179 § 2 欧拉(Euler)方法………………………………………………………………180 ﹒iv﹒ § 3 改进的(Euler)方法……………………………………………………………181 § 4 龙格—库塔(Runge—Kutta)方法…………………………………………… 182 § 5 线性多步法……………………………………………………………………184 § 6 一阶微分方程组与高阶微分方程的数值解法………………………………185 § 7 Matlab 解法………………………………………………………………… 186 习题十五……………………………………………………………………………………191 第十六章 差分方程模型…………………………………………………………………192 § 1 差分方程……………………………………………………………………… 192 § 2 蛛网模型……………………………………………………………………… 195 § 3 商品销售量预测……………………………………………………………… 198 § 4 遗传模型……………………………………………………………………… 200 习题十六……………………………………………………………………………………206 第十七章 马氏链模型……………………………………………………………………207 § 1 随机过程的概念……………………………………………………………… 207 § 2 马尔可夫链…………………………………………………………………… 207 § 3 马尔可夫链的应用…………………………………………………………… 216 习题十七……………………………………………………………………………………216 第十八章 动态优化模型…………………………………………………………………218 § 1 变分法简介…………………………………………………………………… 218 § 2 生产设备的最大经济效益…………………………………………………… 226 习题十八……………………………………………………………………………………229 第十九章 神经网络模型…………………………………………………………………230 § 1 神经网络简介………………………………………………………………… 230 § 2 蠓虫分类问题与多层前馈网络……………………………………………… 232 § 3 处理蠓虫分类的另一种网络方法…………………………………………… 236 习题十九……………………………………………………………………………………238 第二十章 偏微分方程的数值解…………………………………………………………240 § 1 偏微分方程的定解问题……………………………………………………… 240 § 2 偏微分方程的差分解法……………………………………………………… 242 § 3 Matlab 解法……………………………………………………………………247 习题二十……………………………………………………………………………………251 第二十一章 目标规划……………………………………………………………………253 § 1 目标规划的数学模型………………………………………………………… 253 § 2 多目标规划的 Matlab 解法…………………………………………………… 255 习题二十一…………………………………………………………………………………256 第二十二章 模糊数学模型………………………………………………………………257 § 1 模糊数学基本知识…………………………………………………………… 257 § 2 模糊分类问题………………………………………………………………… 263 § 3 最佳方案的模糊决策………………………………………………………… 268 ﹒v﹒ 第二十三章 现代优化算法简介…………………………………………………………271 § 1 现代优化算法简介…………………………………………………………… 271 § 2 模拟退火算法………………………………………………………………… 271 § 3 蚁群算法……………………………………………………………………… 277 第二十四章 时间序列模型………………………………………………………………280 § 1 确定性时间序列分析方法概述……………………………………………… 280 § 2 平稳时间序列模型…………………………………………………………… 284 § 3 ARMA 模型的特性…………………………………………………………… 285 § 4 时间序列建模的基本步骤…………………………………………………… 288 附录一 Matlab入门………………………………………………………………………291 附录二 Matlab在线性代数中的应用……………………………………………………310 附录三 运筹学的 Lingo 软件……………………………………………………………314 附录四 判别分析…………………………………………………………………………317 参考文献……………………………………………………………………………………319 -1- 前 言 今天,人类社会正处在由工业化社会向信息化社会过渡的变革。以数字化为特征的 信息社会有两个显著特点:计算机技术的迅速发展与广泛应用;数学的应用向一切领域 渗透。随着计算机技术的飞速发展,科学计算的作用越来越引起人们的广泛重视,它已 经与科学理论和科学实验并列成为人们探索和研究自然界、人类社会的三大基本方法。 为了适应这种社会的变革,培养和造就出一批又一批适应高度信息化社会具有创新能力 的高素质的工程技术和管理人才,在各高校开设“数学建模”课程,培养学生的科学计 算能力和创新能力,就成为这种新形势下的历史必然。 数学建模是对现实世界的特定对象,为了特定的目的,根据特有的内在规律,对其 进行必要的抽象、归纳、假设和简化,运用适当的数学工具建立的一个数学结构。数学 建模就是运用数学的思想方法、数学的语言去近似地刻画一个实际研究对象,构建一座 沟通现实世界与数学世界的桥梁,并以计算机为工具应用现代计算技术达到解决各种实 际问题的目的。建立一个数学模型的全过程称为数学建模。因此“数学建模”(或数学 实验)课程教学对于开发学生的创新意识,提升人的数学素养,培养学生创造性地应用 数学工具解决实际问题的能力,有着独特的功能。 数学建模过程就是一个创造性的工作过程。人的创新能力首先是创造性思维和具备 创新的思想方法。数学本身是一门理性思维科学,数学教学正是通过各个教学环节对学 生进行严格的科学思维方法的训练,从而引发人的灵感思维,达到培养学生的创造性思 维的能力。同时数学又是一门实用科学,它具有能直接用于生产和实践,解决工程实际 中提出的问题,推动生产力的发展和科学技术的进步。学生参加数学建模活动,首先就 要了解问题的实际背景,深入到具体学科领域的前沿,这就需要学生具有能迅速查阅大 量科学资料,准确获得自己所需信息的能力;同时,不但要求学生必需了解现代数学各 门学科知识和各种数学方法,把所掌握的数学工具创造性地应用于具体的实际问题,构 建其数学结构,还要求学生熟悉各种数学软件,熟练地把现代计算机技术应用于解决当 -2- 前实际问题综合能力,最后还要具有把自己的实践过程和结果叙述成文字的写作能力。 通过数学建模全过程的各个环节,学生们进行着创造性的思维活动,模拟了现代科学研 究过程。通过“数学建模”课程的教学和数学建模活动极大地开发了学生的创造性思维 的能力,培养学生在面对错综复杂的实际问题时,具有敏锐的观察力和洞察力,以及丰 富的想象力。因此,“数学建模”课程在培养学生的创新能力方面有着其它课程不可替 代的作用。 几年的“数学建模”教学实践告诉我们,进行数学建模教学,为学生提供一本内容 丰富,既理论完整又实用的“数学建模”教材,使学生少走弯路尤为重要。这也是我们 编写这本教材的初衷。本教材可以说既是我们多年教学经验的总结,也是我们心血的结 晶。本教材的特点是尽量为学生提供常用的数学方法,并将相应的 Matlab 程序提供给 学生,使学生在进行书中提供的案例的学习中,在自己动手构建数学模型的同时上机进 行数学实验,从而为学生提供数学建模全过程的训练,以便能够达到举一反三,取得事 半功倍的教学效果。司守奎同志编写了全部的 Matlab 程序,参加本书编写的还有毛凯 同志。 全书共二十四章,各章有一定的独立性,这样便于教师和学生按需要进行选择。完 成本教材的教学大约需要 60 学时,其中方法教学与上机实践的比例一般不应少于 1:1 。 一本好的教材需要经过多年的教学实践,反复锤炼。由于我们的经验和时间所限, 书中的错误和纰漏在所难免,敬请同行不吝指正。 编者 2003年 12 月 -1- 第一章 线性规划 §1 线性规划 在人们的生产实践中,经常会遇到如何利用现有资源来安排生产,以取得最大经济 效益的问题。此类问题构成了运筹学的一个重要分支—数学规划,而线性规划(Linear Programming 简记 LP)则是数学规划的一个重要分支。自从 1947 年 G. B. Dantzig 提出 求解线性规划的单纯形方法以来,线性规划在理论上趋向成熟,在实用中日益广泛与深 入。特别是在计算机能处理成千上万个约束条件和决策变量的线性规划问题之后,线性 规划的适用领域更为广泛了,已成为现代管理中经常采用的基本方法之一。 1.1 线性规划的实例与定义 例 1 某机床厂生产甲、乙两种机床,每台销售后的利润分别为 4000 元与 3000 元。 生产甲机床需用 BA、 机器加工,加工时间分别为每台 2 小时和 1 小时;生产乙机床 需用 CBA 、、 三种机器加工,加工时间为每台各一小时。若每天可用于加工的机器时 数分别为 A 机器 10 小时、B 机器 8 小时和C 机器 7 小时,问该厂应生产甲、乙机床各 几台,才能使总利润最大? 上述问题的数学模型:设该厂生产 1x 台甲机床和 2x 乙机床时总利润最大,则 21, xx 应满足 (目标函数) 21 34max xxz += (1) s.t.(约束条件) ⎪ ⎪ ⎩ ⎪⎪ ⎨ ⎧ ≥ ≤ ≤+ ≤+ 0, 7 8 102 21 2 21 21 xx x xx xx (2) 这里变量 21, xx 称之为决策变量,(1)式被称为问题的目标函数,(2)中的几个不等式 是问题的约束条件,记为 s.t.(即 subject to)。由于上面的目标函数及约束条件均为线性 函数,故被称为线性规划问题。 总之,线性规划问题是在一组线性约束条件的限制下,求一线性目标函数最大或最 小的问题。 在解决实际问题时,把问题归结成一个线性规划数学模型是很重要的一步,但往往 也是困难的一步,模型建立得是否恰当,直接影响到求解。而选适当的决策变量,是我 们建立有效模型的关键之一。 1.2 线性规划的 Matlab 标准形式 线性规划的目标函数可以是求最大值,也可以是求最小值,约束条件的不等号可以 是小于号也可以是大于号。为了避免这种形式多样性带来的不便,Matlab 中规定线性 规划的标准形式为 xc x T min s.t. ⎪⎩ ⎪⎨ ⎧ ≤≤ =⋅ ≤ ubxlb beqxAeq bAx 其中 c 和 x 为 n 维列向量, A 、 Aeq 为适当维数的矩阵,b 、beq 为适当维数的列向 量。 -2- 例如线性规划 bAxxc x T ≥ s.t. max 的 Matlab 标准型为 bAxxc x T −≤−− s.t. min 1.3 线性规划问题的解的概念 一般线性规划问题的(数学)标准型为 ∑ = = n j jj xcz 1 max (3) s.t. ⎪⎩ ⎪⎨ ⎧ =≥ ==∑ = njx mibxa j n j ijij ,,2,10 ,,2,1 1 L L (4) 可行解 满足约束条件(4)的 解 ),,,( 21 nxxxx L= ,称为线性规划问题的可行解, 而使目标函数(3)达到最大值的可行解叫最优解。 可行域 所有可行解构成的集合称为问题的可行域,记为 R 。 1.4 线性规划的图解法 0 2 4 6 8 10 0 1 2 3 4 5 6 7 8 9 10 x2 =7 2x1+x2=10 x1 +x2 =8 z=1 2 (2 ,6 ) 图 1 线性规划的图解示意图 图解法简单直观,有助于了解线性规划问题求解的基本原理。我们先应用图解法来 求解例 1。对于每一固定的值 z ,使目标函数值等于 z 的点构成的直线称为目标函数等 位线,当 z 变动时,我们得到一族平行直线。对于例 1,显然等位线越趋于右上方,其 上的点具有越大的目标函数值。不难看出,本例的最优解为 Tx )6,2(* = ,最优目标值 26* =z 。 从上面的图解过程可以看出并不难证明以下断言: (1)可行域 R 可能会出现多种情况。R 可能是空集也可能是非空集合,当 R 非空 时,它必定是若干个半平面的交集(除非遇到空间维数的退化)。R 既可能是有界区域, 也可能是无界区域。 (2)在 R 非空时,线性规划既可以存在有限最优解,也可以不存在有限最优解(其 目标函数值无界)。 -3- (3)若线性规划存在有限最优解,则必可找到具有最优目标函数值的可行域 R 的 “顶点”。 上述论断可以推广到一般的线性规划问题,区别只在于空间的维数。在一般的 n 维 空间中,满足一线性等式 ∑ = = n i ii bxa 1 的点集被称为一个超平面,而满足一线性不等式 ∑ = ≤ n i ii bxa 1 (或 ∑ = ≥ n i ii bxa 1 )的点集被称为一个半空间(其中 ),,( 1 naa L 为一 n 维行 向量, b 为一实数)。若干个半空间的交集被称为多胞形,有界的多胞形又被称为多面 体。易见,线性规划的可行域必为多胞形(为统一起见,空集Φ 也被视为多胞形)。 在一般 n 维空间中,要直接得出多胞形“顶点”概念还有一些困难。二维空间中的顶点 可以看成为边界直线的交点,但这一几何概念的推广在一般 n 维空间中的几何意义并不 十分直观。为此,我们将采用另一途径来定义它。 定义 1 称 n 维空间中的区域 R 为一凸集,若 Rxx ∈∀ 21, 及 )1,0(∈∀λ ,有 Rxx ∈−+ 21 )1( λλ 。 定义 2 设 R 为 n 维空间中的一个凸集, R 中的点 x 被称为 R 的一个极点,若不 存在 Rxx ∈21、 及 )1,0(∈λ ,使得 21 )1( xxx λλ −+= 。 定义 1 说明凸集中任意两点的连线必在此凸集中;而定义 2 说明,若 x 是凸集 R 的一个极点,则 x 不能位于 R 中任意两点的连线上。不难证明,多胞形必为凸集。同 样也不难证明,二维空间中可行域 R 的顶点均为 R 的极点( R 也没有其它的极点)。 1.5 求解线性规划的 Matlab 解法 单纯形法是求解线性规划问题的最常用、最有效的算法之一。这里我们就不介绍 单纯形法,有兴趣的读者可以参看其它线性规划书籍。下面我们介绍线性规划的 Matlab 解法。 Matlab 中线性规划的标准型为 xc x T min s.t. ⎪⎩ ⎪⎨ ⎧ ≤≤ =⋅ ≤ ubxlb beqxAeq bAx 基本函数形式为 linprog(c,A,b),它的返回值是向量 x 的值。还有其它的一些函数调用形 式(在 Matlab 指令窗运行 help linprog 可以看到所有的函数调用形式),如: [x,fval]=linprog(c,A,b,Aeq,beq,LB,UB,X0,OPTIONS) 这里 fval 返回目标函数的值,LB 和 UB 分别是变量 x 的下界和上界, 0x 是 x 的初始值, OPTIONS 是控制参数。 例 2 求解下列线性规划问题 321 532max xxxz −+= s.t. 7321 =++ xxx 1052 321 ≥+− xxx 123 321 ≤++ xxx 0,, 321 ≥xxx -4- 解 (i)编写 M 文件 c=[2;3;-5]; a=[-2,5,-1;1,3,1]; b=[-10;12]; aeq=[1,1,1]; beq=7; x=linprog(-c,a,b,aeq,beq,zeros(3,1)) value=c'*x (ii)将M文件存盘,并命名为example1.m。 (iii)在Matlab指令窗运行example1即可得所求结果。 例3 求解线性规划问题 321 32 min xxxz ++= ⎪⎩ ⎪⎨ ⎧ ≥ ≥+ ≥++ 0,, 623 824 321 21 321 xxx xx xxx 解 编写Matlab程序如下: c=[2;3;1]; a=[1,4,2;3,2,0]; b=[8;6]; [x,y]=linprog(c,-a,-b,[],[],zeros(3,1)) 1.6 可以转化为线性规划的问题 很多看起来不是线性规划的问题也可以通过变换变成线性规划的问题来解决。如: 例4 规划问题为 bAx xxx n ≤ +++ t.s. ||||||min 21 L 其中 T nxxx ][ 1 L= , A 和b 为相应维数的矩阵和向量。 要把上面的问题变换成线性规划问题,只要注意到事实:对任意的 ix ,存在 0, >ii vu 满足 iii vux −= , iii vux +=|| 事实上,我们只要取 2 || ii i xxu += , 2 || ii i xxv −= 就可以满足上面的条件。 这样,记 T nuuu ][ 1 L= , T nvvv ][ 1 L= ,从而我们可以把上面的问题 变成 ∑ = + n i ii vu 1 )(min ⎩ ⎨ ⎧ ≥ ≤− 0, )( t.s. vu bvuA 例 5 |}|max{min iyx ii ε 其中 iii yx −=ε 。 对于这个问题,如果我们取 ||max0 iyi x ε= ,这样,上面的问题就变换成 -5- 0min x 0011 ,, t.s. xyxxyx nn ≤−≤− L 此即我们通常的线性规划问题。 §2 运输问题(产销平衡) 例 6 某商品有 m 个产地、 n 个销地,各产地的产量分别为 maa ,,1 L ,各销地的 需求量分别为 nbb ,,1 L 。若该商品由i 产地运到 j 销地的单位运价为 ijc ,问应该如何调 运才能使总运费最省? 解:引入变量 ijx ,其取值为由i 产地运往 j 销地的该商品数量,数学模型为 ∑∑ == m i n j ijij xc 11 min s.t. ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ ≥ == == ∑ ∑ = = 0 ,,2,1, ,,1, 1 1 ij m i jij n j iij x njbx miax L L 显然是一个线性规划问题,当然可以用单纯形法求解。 对产销平衡的运输问题,由于有以下关系式存在: ∑∑∑∑∑∑ ====== =⎟ ⎠ ⎞⎜ ⎝ ⎛=⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛= m i i n j n j m i ij m i n j ijj axxb 111111 其约束条件的系数矩阵相当特殊,可用比较简单的计算方法,习惯上称为表上作业法(由 康托洛维奇和希奇柯克两人独立地提出,简称康—希表上作业法)。 §3 指派问题 3.1 指派问题的数学模型 例 7 拟分配 n 人去干 n 项工作,每人干且仅干一项工作,若分配第i 人去干第 j 项工作,需花费 ijc 单位时间,问应如何分配工作才能使工人花费的总时间最少? 容易看出,要给出一个指派问题的实例,只需给出矩阵 )( ijcC = ,C 被称为指派 问题的系数矩阵。 引入变量 ijx ,若分配i 干 j 工作,则取 1=ijx ,否则取 0=ijx 。上述指派问题的 数学模型为 ∑∑ == n i n j ijij xc 11 min s.t. ∑ = = n j ijx 1 1 -6- ∑ = = n i ijx 1 1 1 0或=ijx 上述指派问题的可行解可以用一个矩阵表示,其每行每列均有且只有一个元素为 1,其余元素均为 0;可以用 n,,1 L 中的一个置换表示。 问题中的变量只能取 0 或 1,从而是一个 0-1 规划问题。一般的 0-1 规划问题求解 极为困难。但指派问题并不难解,其约束方程组的系数矩阵十分特殊(被称为全单位模 矩阵,其各阶非零子式均为 1± ),其非负可行解的分量只能取 0 或 1,故约束 10或=ijx 可改写为 0≥ijx 而不改变其解。此时,指派问题被转化为一个特殊的运输问题,其中 nm = , 1== ji ba 。 3.2 求解指派问题的匈牙利算法 由于指派问题的特殊性,又存在着由匈牙利数学家 Konig 提出的更为简便的解法 —匈牙利算法。算法主要依据以下事实:如果系数矩阵 )( ijcC = 一行(或一列)中每 一元素都加上或减去同一个数,得到一个新矩阵 )( ijbB = ,则以C 或 B 为系数矩阵的 指派问题具有相同的最优指派。 例 8 求解指派问题,其系数矩阵为 ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = 16221917 17182224 18192117 22191516 C 解 将第一行元素减去此行中的最小元素 15,同样,第二行元素减去 17,第三行 元素减去 17,最后一行的元素减去 16,得 ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = 0631 0157 1240 7401 1B 再将第 3 列元素各减去 1,得 ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = * * * * 2 0531 0057 1140 7301 B 以 2B 为系数矩阵的指派问题有最优指派 ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ 4312 4321 由等价性,它也是例 7 的最优指派。 有时问题会稍复杂一些。 例 9 求解系数矩阵C 的指派问题 -7- ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = 6107104 10661415 121412177 66698 979712 C 解:先作等价变换如下 ∨ ∨ ∨ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ → ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − − − − − 26360 40*089 57510*0 0*0032 202*05 6107104 10661415 121412177 66698 979712 4 6 7 6 7 容易看出,从变换后的矩阵中只能选出四个位于不同行不同列的零元素,但 5=n , 最优指派还无法看出。此时等价变换还可进行下去。步骤如下: (1) 对未选出 0 元素的行打 ∨ ; (2) 对 ∨ 行中 0 元素所在列打 ∨ ; (3) 对 ∨ 列中选中的 0 元素所在行打 ∨ ; 重复(2)、(3)直到无法再打 ∨ 为止。 可以证明,若用直线划没有打 ∨ 的行与打 ∨ 的列,就得到了能够覆盖住矩阵中所 有零元素的最少条数的直线集合,找出未覆盖的元素中的最小者,令 ∨ 行元素减去此 数, ∨ 列元素加上此数,则原先选中的 0 元素不变,而未覆盖元素中至少有一个已转 变为 0,且新矩阵的指派问题与原问题也等价。上述过程可反复采用,直到能选取出足 够的 0 元素为止。例如,对例 5 变换后的矩阵再变换,第三行、第五行元素减去 2,第 一列元素加上 2,得 ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ 04140 400811 35380 00034 20207 现在已可看出,最优指派为 ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ 53142 54321 。 §4 对偶理论与灵敏度分析 4.1 原始问题和对偶问题 考虑下列一对线性规划模型: xcTmax s.t. 0, ≥≤ xbAx (P) 和 ybTmin s.t. 0, ≥≥ ycyAT (D) -8- 称(P)为原始问题,(D)为它的对偶问题。 不太严谨地说,对偶问题可被看作是原始问题的“行列转置”: (1) 原始问题中的第 j 列系数与其对偶问题中的第 j 行的系数相同; (2) 原始目标函数的各个系数行与其对偶问题右侧的各常数列相同; (3) 原始问题右侧的各常数列与其对偶目标函数的各个系数行相同; (4) 在这一对问题中,不等式方向和优化方向相反。 考虑线性规划: 0,s.t.min ≥= xbAxxcT 把其中的等式约束变成不等式约束,可得 0, s.t. min ≥⎥⎦ ⎤ ⎢⎣ ⎡ −≥⎥⎦ ⎤ ⎢⎣ ⎡ − xb bxA AxcT 它的对偶问题是 [] [ ]cy yAAy ybb TTTT ≤⎥⎦ ⎤ ⎢⎣ ⎡−⎥⎦ ⎤ ⎢⎣ ⎡− 2 1 2 1 s.t.max 其中 1y 和 2y 分别表示对应于约束 bAx ≥ 和 bAx −≥− 的对偶变量组。令 21 yyy −= , 则上式又可写成 cyAyb TT ≤ s.t. max 原问题和对偶的对偶约束之间的关系: min max ⎪⎩ ⎪⎨ ⎧ ≤ ≥ 无限制 变量 0 0 ⎪⎩ ⎪⎨ ⎧ = ≥ ≤ 0 0 0 行约束 ⎪⎩ ⎪⎨ ⎧ = ≤ ≥ 0 0 0 行约束 ⎪⎩ ⎪⎨ ⎧ ≤ ≥ 无限制 变量 0 0 4.2 对偶问题的基本性质 1o 对称性:对偶问题的对偶是原问题。 2o 弱对偶性:若 x 是原问题的可行解, y 是对偶问题的可行解。则存在 ybxc TT ≤ 。 3o 无界性:若原问题(对偶问题)为无界解,则其对偶问题(原问题)无可行解。 4o 可行解是最优解时的性质:设 xˆ 是原问题的可行解, yˆ 是对偶问题的可行解, 当 ybxc TT ˆˆ = 时, yx ˆ,ˆ 是最优解。 5o 对偶定理:若原问题有最优解,那么对偶问题也有最优解;且目标函数值相同。 6o 互补松弛性:若 yx ˆ,ˆ 分别是原问题和对偶问题的最优解,则 0)ˆ(ˆ,0)ˆ(ˆ =−=− cyAxbxAy TTT 例 10 已知线性规划问题 54321 32532min xxxxx ++++=ω s.t. 432 54321 ≥++++ xxxxx -9- 332 54321 ≥+++− xxxxx 5,,2,1,0 L=≥ jx j 已知其对偶问题的最优解为 5;5 3,5 4 * 2 * 1 === zyy 。试用对偶理论找出原问题的最优 解。 解 先写出它的对偶问题 21 34max yyz += 22 21 ≤+ yy ① 321 ≤− yy ② 532 31 ≤+ yy ③ 221 ≤+ yy ④ 33 21 ≤+ yy ⑤ 0, 21 ≥yy 将 * 2 * 1 , yy 的值代入约束条件,得②,③,④为严格不等式;由互补松弛性得 0* 4 * 3 * 2 === xxx 。因 0, * 2 * 1 >yy ;原问题的两个约束条件应取等式,故有 43 * 5 * 1 =+ xx 32 * 5 * 1 =+ xx 求解后得到 1,1 * 5 * 1 == xx ;故原问题的最优解为 5;]'10001[ ** == ωX 。 4.3 灵敏度分析 在以前讨论线性规划问题时,假定 jiij cba ,, 都是常数。但实际上这些系数往往是估 计值和预测值。如市场条件一变, jc 值就会变化; ija 往往是因工艺条件的改变而改变; ib 是根据资源投入后的经济效果决定的一种决策选择。因此提出这样两个问题:当这 些系数有一个或几个发生变化时,已求得的线性规划问题的最优解会有什么变化;或者 这些系数在什么范围内变化时,线性规划问题的最优解或最优基不变。这里我们就不讨 论了。 4.4 参数线性规划 参数线性规划是研究 jiij cba ,, 这些参数中某一参数连续变化时,使最优解发生变化 的各临界点的值。即把某一参数作为参变量,而目标函数在某区间内是这参变量的线性 函数,含这参变量的约束条件是线性等式或不等式。因此仍可用单纯形法和对偶单纯形 法进行分析参数线性规划问题。 §5 投资的收益和风险 5.1 问题提出 市场上有 n 种资产 is ( ni ,,2,1 L= )可以选择,现用数额为 M 的相当大的资金 作一个时期的投资。这 n 种资产在这一时期内购买 is 的平均收益率为 ir ,风险损失率为 iq ,投资越分散,总的风险越少,总体风险可用投资的 is 中最大的一个风险来度量。 -10- 购买 is 时要付交易费,(费率 ip ),当购买额不超过给定值 iu 时,交易费按购买 iu 计算。另外,假定同期银行存款利率是 0r ,既无交易费又无风险。( %50 =r ) 已知 4=n 时相关数据如表 1。 表 1 is ir (%) iq ip (%) iu (元) 1s 28 2.5 1 103 2s 21 1.5 2 198 3s 23 5.5 4.5 52 4s 25 2.6 6.5 40 试给该公司设计一种投资组合方案,即用给定资金 M ,有选择地购买若干种资产 或存银行生息,使净收益尽可能大,使总体风险尽可能小。 5.2 符号规定和基本假设 符号规定: is :第i 种投资项目,如股票,债券 iii qpr ,, :分别为 is 的平均收益率,交易费率,风险损失率 iu : is 的交易定额 0r :同期银行利率 ix :投资项目 is 的资金 a :投资风险度 Q :总体收益 基本假设: 1. 投资数额 M 相当大,为了便于计算,假设 1=M ; 2. 投资越分散,总的风险越小; 3. 总体风险用投资项目 is 中最大的一个风险来度量; 4. n 种资产 is 之间是相互独立的; 5. 在投资的这一时期内, iii qpr ,, , 0r 为定值,不受意外因素影响; 6. 净收益和总体风险只受 iii qpr ,, 影响,不受其它因素干扰。 5.3 模型的分析与建立 1. 总体风险用所投资的 is 中最大的一个风险来衡量,即 },,2,1|max{ nixq ii L= 2.购买 is 所付交易费是一个分段函数,即 交易费 ⎩ ⎨ ⎧ ≤ >= iiii iiii uxup uxxp , , 而题目所给定的定值 iu (单位:元)相对总投资 M 很少, iiup 更小,可以忽略不 计,这样购买 is 的净收益为 iii xpr )( − 。 3. 要使净收益尽可能大,总体风险尽可能小,这是一个多目标规划模型: -11- 目标函数为 ⎪⎩ ⎪⎨ ⎧ −∑ = }max{min )(max 0 ii n i iii xq xpr 约束条件为 ⎪⎩ ⎪⎨ ⎧ =≥ =+∑ = nix Mxp i n i ii ,,1,0,0 )1( 0 L 4. 模型简化 a) 在实际投资中,投资者承受风险的程度不一样,若给定风险一个界限 a ,使最 大的一个风险 aM xq ii ≤ ,可找到相应的投资方案。这样把多目标规划变成一个目标的线 性规划。 模型一 固定风险水平,优化收益 ∑ = − n i iii xpr 0 )(max s.t. ⎪ ⎪ ⎩ ⎪⎪ ⎨ ⎧ =≥=+ ≤ ∑ = n i iii ii nixMxp aM xq 0 ,,1,0,0,)1( L b) 若投资者希望总盈利至少达到水平 k 以上,在风险最小的情况下寻求相应的投 资组合。 模型二 固定盈利水平,极小化风险 }}{{maxmin ii xq s.t. ⎪ ⎪ ⎩ ⎪⎪ ⎨ ⎧ =≥=+ ≥− ∑ ∑ = = n i iii n i iii nixMxp kxpr 0 0 ,,1,0,0,)1( )( L c) 投资者在权衡资产风险和预期收益两方面时,希望选择一个令自己满意的投资 组合。因此对风险、收益分别赋予权重 s( 10 ≤< s )和 )1( s− ,s 称为投资偏好系数。 模型三 ∑ = −−− n i iiiii xprsxqs 0 )()1(}}{max{min s.t. ∑ = =≥=+ n i iii nixMxp 0 ,,2,1,0,0,)1( L 5.4 模型一的求解 模型一为: Txxxxxf ),,,,)(185.0,185.0,19.0,27.0,05.0(min 43210−−−−−= -12- s.t. ⎪⎪ ⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎪ ⎨ ⎧ =≥ ≤ ≤ ≤ ≤ =++++ )4,,1,0(0 026.0 055.0 015.0 025.0 1065.1045.102.101.1 4 3 2 1 43210 Lix ax ax ax ax xxxxx i 由于 a 是任意给定的风险度,到底怎样没有一个准则,不同的投资者有不同的风险 度。我们从 0=a 开始,以步长 001.0=Δa 进行循环搜索,编制程序如下: clc,clear a=0; hold on while a<0.05 c=[-0.05,-0.27,-0.19,-0.185,-0.185]; A=[zeros(4,1),diag([0.025,0.015,0.055,0.026])]; b=a*ones(4,1); Aeq=[1,1.01,1.02,1.045,1.065]; beq=1; LB=zeros(5,1); [x,Q]=linprog(c,A,b,Aeq,beq,LB); Q=-Q; plot(a,Q,'*r'); a=a+0.001; end xlabel('a'),ylabel('Q') 5.5 结果分析 1. 风险大,收益也大。 2.当投资越分散时,投资者承担的风险越小,这与题意一致。即: 冒险的投资者会出现集中投资的情况,保守的投资者则尽量分散投资。 3.在 006.0=a 附近有一个转折点,在这一点左边,风险增加很少时,利润增长 很快。在这一点右边,风险增加很大时,利润增长很缓慢,所以对于风险和收益没有特 殊偏好的投资者来说,应该选择曲线的拐点作为最优投资组合,大约是 %6.0=a , %20=Q ,所对应投资方案为: 风险度 006.0=a ,收益 2019.0=Q , 00 =x , 24.01 =x , 4.02 =x , 1091.03 =x , 2212.04 =x 。 习 题 一 1.试将下述问题改写成线性规划问题: ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ⎟ ⎠ ⎞⎜ ⎝ ⎛ ∑∑∑ === m i iin m i ii m i iix xaxaxa i 11 2 1 1 ,,,minmax L ⎩ ⎨ ⎧ =≥ =+++ mix xxx i m ,,1,0 1st. 21 L L 2.试将下列问题改写成线性规划问题: -13- ∑ = = n j jj xcz 1 ||max ⎪⎩ ⎪⎨ ⎧ ==∑ = 取值无约束j n j ijij x mibxa 1 ),,2,1( st. L 3.线性回归是一种常用的数理统计方法,这个方法要求对图上的一系列点 ),(,),,(),,( 2211 nn yxyxyx L 选配一条合适的直线拟合。方法通常是先定直线方程为 bxay += ,然后按某种准则求定 ba, 。通常这个准则为最小二乘法,但也可用其他准 则。试根据以下准则建立这个问题的线性规划模型: ∑ = +− n i ii bxay 1 |)(|min 4.某厂生产三种产品 I,II,III。每种产品要经过 BA, 两道工序加工。设该厂有 两种规格的设备能完成 A 工序,它们以 21, AA 表示;有三种规格的设备能完成 B 工序, 它们以 321 ,, BBB 表示。产品 I 可在 BA, 任何一种规格设备上加工。产品 II 可在任何规 格的 A 设备上加工,但完成 B 工序时,只能在 1B 设备上加工;产品 III 只能在 2A 与 2B 设备上加工。已知在各种机床设备的单件工时,原材料费,产品销售价格,各种设备有 效台时以及满负荷操作时机床设备的费用如表 2,求安排最优的生产计划,使该厂利润 最大。 表 2 产 品 设备 I II III 设备有效台时 满负荷时的 设备费用(元) 1A 2A 1B 2B 3B 5 7 6 4 7 10 9 8 12 11 6000 10000 4000 7000 4000 300 321 250 783 200 原料费(元/件) 单 价(元/件) 0.25 1.25 0.35 2.00 0.50 2.80 5.有四个工人,要指派他们分别完成 4 项工作,每人做各项工作所消耗的时间如 表 3。 表 3 工作 工人 A B C D 甲 15 18 21 24 乙 19 23 22 18 丙 26 17 16 19 丁 19 21 23 17 问指派哪个人去完成哪项工作,可使总的消耗时间为最小? -14- 6.某战略轰炸机群奉命摧毁敌人军事目标。已知该目标有四个要害部位,只要摧 毁其中之一即可达到目的。为完成此项任务的汽油消耗量限制为 48000 升、重型炸弹 48 枚、轻型炸弹 32 枚。飞机携带重型炸弹时每升汽油可飞行 2 千米,带轻型炸弹时每 升汽油可飞行 3 千米。又知每架飞机每次只能装载一枚炸弹,每出发轰炸一次除来回路 程汽油消耗(空载时每升汽油可飞行 4 千米)外,起飞和降落每次各消耗 100 升。有关 数据如表 4 所示。 表 4 摧毁可能性 要害部位 离机场距离 (千米) 每枚重型弹 每枚轻型弹 1 2 3 4 450 480 540 600 0.10 0.20 0.15 0.25 0.08 0.16 0.12 0.20 为了使摧毁敌方军事目标的可能性最大,应如何确定飞机轰炸的方案,要求建立这个问 题的线性规划模型。 7.用 Matlab 求解下列线性规划问题: 3213max xxxz −−= s.t. ⎪ ⎪ ⎩ ⎪⎪ ⎨ ⎧ ≥ =+− ≥++− ≤+− 0,, 12 324 112 321 31 321 321 xxx xx xxx xxx 8.用 Matlab 求解下列规划问题: ||4||3||2||min 4321 xxxxz +++= s.t. 04321 =+−− xxxx 13 4321 =−+− xxxx 2 132 4321 −=+−− xxxx 9.一架货机有三个货舱:前舱、中仓和后舱。三个货舱所能装载的货物的最大重 量和体积有限制如表 5 所示。并且为了飞机的平衡,三个货舱装载的货物重量必须与其 最大的容许量成比例。 表 5 货舱数据 前舱 中仓 后舱 重量限制(吨) 10 16 8 体积限制(立方米) 6800 8700 5300 现有四类货物用该货机进行装运,货物的规格以及装运后获得的利润如表 6。 表 6 货物规格及利润表 重量(吨) 空间(立方米/吨) 利润(元/吨) 货物 1 18 480 3100 货物 2 15 650 3800 货物 3 23 580 3500 货物 4 12 390 2850 -15- 假设: (1)每种货物可以无限细分; (2)每种货物可以分布在一个或者多个货舱内; (3)不同的货物可以放在同一个货舱内,并且可以保证不留空隙。 问应如何装运,使货机飞行利润最大? -16- 第二章 整数规划 §1 概论 1.1 定义 规划中的变量(部分或全部)限制为整数时,称为整数规划。若在线性规划模型中, 变量限制为整数,则称为整数线性规划。目前所流行的求解整数规划的方法,往往只适 用于整数线性规划。目前还没有一种方法能有效地求解一切整数规划。 1.2 整数规划的分类 如不加特殊说明,一般指整数线性规划。对于整数线性规划模型大致可分为两类: 1o 变量全限制为整数时,称纯(完全)整数规划。 2o 变量部分限制为整数的,称混合整数规划。 1.2 整数规划特点 (i) 原线性规划有最优解,当自变量限制为整数后,其整数规划解出现下述情况: ①原线性规划最优解全是整数,则整数规划最优解与线性规划最优解一致。 ②整数规划无可行解。 例 1 原线性规划为 21min xxz += 0,0,542 2121 ≥≥=+ xxxx 其最优实数解为: 4 5min,4 5,0 21 === zxx 。 ③有可行解(当然就存在最优解),但最优解值变差。 例 2 原线性规划为 21min xxz += 0,0,642 2121 ≥≥=+ xxxx 其最优实数解为: 2 3min,2 3,0 21 === zxx 。 若限制整数得: 2min,1,1 21 === zxx 。 (ii) 整数规划最优解不能按照实数最优解简单取整而获得。 1.3 求解方法分类: (i)分枝定界法—可求纯或混合整数线性规划。 (ii)割平面法—可求纯或混合整数线性规划。 (iii)隐枚举法—求解“0-1”整数规划: ①过滤隐枚举法; ②分枝隐枚举法。 (iv)匈牙利法—解决指派问题(“0-1”规划特殊情形)。 (v)蒙特卡洛法—求解各种类型规划。 下面将简要介绍常用的几种求解整数规划的方法。 §2 分枝定界法 对有约束条件的最优化问题(其可行解为有限数)的所有可行解空间恰当地进行系 统搜索,这就是分枝与定界内容。通常,把全部可行解空间反复地分割为越来越小的子 集,称为分枝;并且对每个子集内的解集计算一个目标下界(对于最小值问题),这称 为定界。在每次分枝后,凡是界限超出已知可行解集目标值的那些子集不再进一步分枝, -17- 这样,许多子集可不予考虑,这称剪枝。这就是分枝定界法的主要思路。 分枝定界法可用于解纯整数或混合的整数规划问题。在本世纪六十年代初由 Land Doig 和 Dakin 等人提出的。由于这方法灵活且便于用计算机求解,所以现在它已是解 整数规划的重要方法。目前已成功地应用于求解生产进度问题、旅行推销员问题、工厂 选址问题、背包问题及分配问题等。 设有最大化的整数规划问题 A ,与它相应的线性规划为问题 B ,从解问题 B 开始, 若其最优解不符合 A 的整数条件,那么 B 的最优目标函数必是 A 的最优目标函数 *z 的 上界,记作 z ;而 A 的任意可行解的目标函数值将是 *z 的一个下界 z 。分枝定界法就 是将 B 的可行域分成子区域的方法。逐步减小 z 和增大 z ,最终求到 *z 。现用下例来 说明: 例 3 求解下述整数规划 21 9040Max xxz += ⎪⎩ ⎪⎨ ⎧ ≥ ≤+ ≤+ 且为整数0, 70207 5679 21 21 21 xx xx xx 解 (i)先不考虑整数限制,即解相应的线性规划 B ,得最优解为: 355.8779,8168.1,8092.4 21 === zxx 可见它不符合整数条件。这时 z 是问题 A 的最优目标函数值 *z 的上界,记作 z 。而 0,0 21 == xx 显然是问题 A 的一个整数可行解,这时 0=z ,是 *z 的一个下界,记作 z , 即 3560 * ≤≤ z 。 (ii)因 为 21, xx 当前均为非整数,故不满足整数要求,任选一个进行分枝。设选 1x 进行分枝,把可行集分成 2 个子集: 44.8092][1 =≤x , 514.8092][1 =+≥x 因为 4 与 5 之间无整数,故这两个子集的整数解必与原可行集合整数解一致。这 一步称为分枝。这两个子集的规划及求解如下: 问题 1B : 21 9040Max xxz += ⎪⎩ ⎪⎨ ⎧ ≥≤≤ ≤+ ≤+ 0,40 70207 5679 21 21 21 xx xx xx 最优解为: 349,1.2,0.4 121 === zxx 。 问题 2B : 21 9040Max xxz += ⎪⎩ ⎪⎨ ⎧ ≥≥ ≤+ ≤+ 0,5 70207 5679 21 21 21 xx xx xx 最优解为: 4.341,57.1,0.5 121 === zxx 。 再定界: 3490 * ≤≤ z 。 (iii)对问题 1B 再进行分枝得问题 11B 和 12B ,它们的最优解为 -18- 340,2,4: 112111 === zxxB 327.14,00.3x1.43,: 122112 === zxB 再定界: 341340 * ≤≤ z ,并将 12B 剪枝。 (iv)对问题 2B 再进行分枝得问题 21B 和 22B ,它们的最优解为 083,00.1x5.44,: 222121 === zxB 22B 无可行解。 将 2221,BB 剪枝。 于是可以断定原问题的最优解为: 340,2,4 * 21 === zxx 从以上解题过程可得用分枝定界法求解整数规划(最大化)问题的步骤为: 开始,将要求解的整数规划问题称为问题 A ,将与它相应的线性规划问题称为问 题 B 。 (i)解问题 B 可能得到以下情况之一: (a) B 没有可行解,这时 A 也没有可行解,则停止. (b) B 有最优解,并符合问题 A 的整数条件, B 的最优解即为 A 的最优解,则 停止。 (c) B 有最优解,但不符合问题 A 的整数条件,记它的目标函数值为 z 。 (ii)用观察法找问题 A 的一个整数可行解,一般可取 njx j ,,1,0 L== ,试探, 求得其目标函数值,并记作 z 。以 *z 表示问题 A 的最优目标函数值;这时有 zzz ≤≤ * 进行迭代。 第一步:分枝,在 B 的最优解中任选一个不符合整数条件的变量 jx ,其值为 jb , 以 ][ jb 表示小于 jb 的最大整数。构造两个约束条件 ][ jj bx ≤ 和 1][ +≥ jj bx 将这两个约束条件,分别加入问题 B ,求两个后继规划问题 1B 和 2B 。不考虑整数条件 求解这两个后继问题。 定界,以每个后继问题为一分枝标明求解的结果,与其它问题的解的结果中,找出 最优目标函数值最大者作为新的上界 z 。从已符合整数条件的各分支中,找出目标函数 值为最大者作为新的下界 z ,若无作用 z 不变。 第二步:比较与剪枝,各分枝的最优目标函数中若有小于 z 者,则剪掉这枝,即 以后不再考虑了。若大于 z ,且不符合整数条件,则重复第一步骤。一直到最后得到 zz =* 为止。得最优整数解 njx j ,,1,* L= 。 §3 10 − 型整数规划 10 − 型整数规划是整数规划中的特殊情形,它的变量 jx 仅取值 0 或 1。这时 jx 称 为 10 − 变量,或称二进制变量。 jx 仅取值 0 或 1 这个条件可由下述约束条件: 10 ≤≤ jx ,整数 -19- 所代替,是和一般整数规划的约束条件形式一致的。在实际问题中,如果引入 10 − 变 量,就可以把有各种情况需要分别讨论的线性规划问题统一在一个问题中讨论了。我们 先介绍引入 10 − 变量的实际问题,再研究解法。 3.1 引入 10 − 变量的实际问题 3.1.1 投资场所的选定——相互排斥的计划 例 4 某公司拟在市东、西、南三区建立门市部。拟议中有 7 个位置(点) )7,,2,1( L=iAi 可供选择。规定 在东区。由 321 ,, AAA 三个点中至多选两个; 在西区。由 54 , AA 两个点中至少选一个; 在南区,由 76 , AA 两个点中至少选一个。 如选用 iA 点,设备投资估计为 ib 元,每年可获利润估计为 ic 元,但投资总额不能 超过 B 元。问应选择哪几个点可使年利润为最大? 解题时先引入 10 − 变量 )7,,2,1( L=ixi 令 ⎩ ⎨ ⎧= .0 ,1 点没被选中当 点被选中当 , , iA iA ix 7,,2,1 L=i . 于是问题可列写成: i i i xcz ∑ = = 7 1 Max ⎪ ⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ =≥+ ≥+ ≤++ ≤∑ = 10,1 1 2 76 54 321 7 1 或i i ii xxx xx xxx Bxb 3.1.2 相互排斥的约束条件 有两个相互排斥的约束条件 2445 21 ≤+ xx 或 4537 21 ≤+ xx 。 为了统一在一个问题中,引入 10 − 变量 y ,则上述约束条件可改写为: ⎪⎩ ⎪⎨ ⎧ = −+≤+ +≤+ 10 )1(4537 2445 21 21 或y Myxx yMxx 其中 M 是充分大的数。 约束条件 01 =x 或 800500 1 ≤≤ x 可改写为 ⎩ ⎨ ⎧ = ≤≤ 10 800500 1 或y yxy -20- 如果有 m 个互相排斥的约束条件: mibxaxa inini ,,2,111 LL =≤++ 为了保证这 m 个约束条件只有一个起作用,我们引入 m 个 10 − 变量 ),,2,1( miyi L= 和一个充分大的常数 M ,而下面这一组 1+m 个约束条件 miMybxaxa iinini ,,2,111 LL =+≤++ (1) 11 −=++ myy mL (2) 就合于上述的要求。这是因为,由于(2), m 个 iy 中只有一个能取 0 值,设 0* =iy , 代入(1),就只有 *ii = 的约束条件起作用,而别的式子都是多余的。 3.1.3 关于固定费用的问题(Fixed Cost Problem) 在讨论线性规划时,有些问题是要求使成本为最小。那时总设固定成本为常数,并 在线性规划的模型中不必明显列出。但有些固定费用(固定成本)的问题不能用一般线 性规划来描述,但可改变为混合整数规划来解决,见下例。 例 5 某工厂为了生产某种产品,有几种不同的生产方式可供选择,如选定的生产 方式投资高(选购自动化程度高的设备),由于产量大,因而分配到每件产品的变动成 本就降低;反之,如选定的生产方式投资低,将来分配到每件产品的变动成本可能增加。 所以必须全面考虑。今设有三种方式可供选择,令 jx 表示采用第 j 种方式时的产量; jc 表示采用第 j 种方式时每件产品的变动成本; jk 表示采用第 j 种方式时的固定成本。 为了说明成本的特点,暂不考虑其它约束条件。采用各种生产方式的总成本分别为 ⎪⎩ ⎪⎨ ⎧ = >+ = 0 ,0 0, j jjjj j x xxck P 当 当 3,2,1=j . 在构成目标函数时,为了统一在一个问题中讨论,现引入 10 − 变量 jy ,令 ⎪⎩ ⎪⎨ ⎧ = = > .0 0 ,0 ,1 时种生产方式,即当不采用第 时,种生产方式,即当采用第 jxj jxj jy (3) 于是目标函数 )()()(min 333322221111 xcykxcykxcykz +++++= (3)式这个规定可表为下述 3 个线性约束条件: 3,2,1, =≤≤ jMyxy jjjε (4) 其中ε 是一个充分小的正常数,M 是个充分大的正常数。(4)式说明,当 0>jx 时 jy 必须为 1;当 0=jx 时只有 jy 为 0 时才有意义,所以(4)式完全可以代替(3)式。 3.2 10 − 型整数规划解法之一(过滤隐枚举法) 解 10 − 型整数规划最容易想到的方法,和一般整数规划的情形一样,就是穷举法, 即检查变量取值为 0 或 1 的每一种组合,比较目标函数值以求得最优解,这就需要检查 变量取值的 n2 个组合。对于变量个数 n 较大(例如 100>n ),这几乎是不可能的。因 此常设计一些方法,只检查变量取值的组合的一部分,就能求到问题的最优解。这样的 方法称为隐枚举法(Implicit Enumeration),分枝定界法也是一种隐枚举法。当然,对 -21- 有些问题隐枚举法并不适用,所以有时穷举法还是必要的。 下面举例说明一种解 10 − 型整数规划的隐枚举法。 例 6 321 523Max xxxz +−= ⎪ ⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ = ≤+ ≤+ ≤++ ≤−+ 10,, 64 3 44 22 321 32 21 321 321 或xxx xx xx xxx xxx 求解思路及改进措施: (i) 先试探性求一个可行解,易看出 )0,0,1(),,( 321 =xxx 满足约束条件,故为一 个可行解,且 3=z 。 (ii) 因为是求极大值问题,故求最优解时,凡是目标值 3x ,这个条件可以表示为 0)500( 21 =− xx (15) 同理,只有当以 8 千元/吨的价格购买 5002 =x (吨)时,才能以 6 千元/吨的价格购买 )0(3 >x ,于是 0)500( 32 =− xx (16) 此外, 321 ,, xxx 的取值范围是 500,,0 321 ≤≤ xxx (17) 由于有非线性约束(15)、(16),因而(7)~(17)构成非线性规划模型。将该模 型输入 LINGO 软件如下: model: sets: var1/1..4/:y; !这里y(1)=x11,y(2)=x21,y(3)=x12,y(4)=x22; var2/1..3/:x,c; endsets max=4.8*(y(1)+y(2))+5.6*(y(3)+y(4))-@sum(var2:c*x); y(1)+y(3)<@sum(var2:x)+500; y(2)+y(4)<1000; 0.5*(y(1)-y(2))>0; 0.4*y(3)-0.6*y(4)>0; (x(1)-500)*x(2)=0; (x(2)-500)*x(3)=0; @for(var2:@bnd(0,x,500)); data: c=10 8 6; enddata end 可以用菜单命令“LINGO|Options”在“Global Solver”选项卡上启动全局优化选 型,并运行上述程序求得全局最有解:购买 1000 吨原油 A ,与库存的 500 吨原油 A 和 1000 吨原油 B 一起,共生产 2500 吨汽油乙,利润为 5000(千元)。 (2)解法二 引入 0-1 变量将(15)和(16)转化为线性约束。 令 11 =z , 12 =z , 13 =z 分别表示以 10 千元/吨、8 千元/吨、6 千元/吨的价格采 购原油 A ,则约束(15)和(16)可以替换为 -26- 112 500500 zxz ≤≤ , (18) 223 500500 zxz ≤≤ , (19) 33 500zx ≤ , (20) 0,, 321 =zzz 或 1 (21) 式(7)~(14),式(17)~(21)构成混合整数线性规划模型,将它输入 LINGO 软 件如下: model: sets: var1/1..4/:y; !这里y(1)=x11,y(2)=x21,y(3)=x12,y(4)=x22; var2/1..3/:x,z,c; endsets max=4.8*(y(1)+y(2))+5.6*(y(3)+y(4))-@sum(var2:c*x); y(1)+y(3)<@sum(var2:x)+500; y(2)+y(4)<1000; 0.5*(y(1)-y(2))>0; 0.4*y(3)-0.6*y(4)>0; @for(var1(i)|i #lt# 3:500*z(i+1)0; 0.4*y(3)-0.6*y(4)>0; w(1)δ ,使 ),0(),()( δ∈∀<+ txftpxf , 称向量 p 是 f 在点 x 处的下降方向。 设 0, ≠∈ pRx n ,若存在 0>t ,使 Ktpx ∈+ , 称向量 p 是点 x 处关于 K 的可行方向。 一个向量 p ,若既是函数 f 在点 x 处的下降方向,又是该点关于区域 K 的可行方 向,称之为函数 f 在点 x 处关于 K 的可行下降方向。 现在,我们给出用基本迭代格式(1)求解(NP)的一般步骤如下: 0° 选取初始点 0x ,令 0:=k 。 1° 构造搜索方向,依照一定规划,构造 f 在点 kx 处关于 K 的可行下降方向作为 搜索方向 kp 。 2° 寻求搜索步长。以 kx 为起点沿搜索方向 kp 寻求适当的步长 kt ,使目标函数值 有某种意义的下降。 3° 求出下一个迭代点。按迭代格式(1)求出 k k kk ptxx +=+1 。 若 1+kx 已满足某种终止条件,停止迭代。 4° 以 1+kx 代替 kx ,回到 1°步。 1.5 凸函数、凸规划 设 )(xf 为定义在 n 维欧氏空间 )(nE 中某个凸集 R 上的函数,若对任何实数 )10( <<αα 以及 R 中的任意两点 )1(x 和 )2(x ,恒有 )()1()())1(( )2()1()2()1( xfxfxxf αααα −+≤−+ 则称 )(xf 为定义在 R 上的凸函数。 若对每一个 )10( << αα 和 Rxx ∈≠ )2()1( 恒有 )()1()())1(( )2()1()2()1( xfxfxxf αααα −+<−+ 则称 )(xf 为定义在 R 上的严格凸函数。 考虑非线性规划 ⎪⎩ ⎪⎨ ⎧ =≤= ∈ },,2,1,0)(|{ )( min ljxgxR xf j Rx L 假定其中 )(xf 为凸函数, ),,2,1)(( ljxg j L= 为凸函数,这样的非线性规划称为 凸规划。 可以证明,凸规划的可行域为凸集,其局部最优解即为全局最优解,而且其最优 解的集合形成一个凸集。当凸规划的目标函数 )(xf 为严格凸函数时,其最优解必定唯 一(假定最优解存在)。由此可见,凸规划是一类比较简单而又具有重要理论意义的非 -36- 线性规划。 §2 无约束问题 2.1 一维搜索方法 当用迭代法求函数的极小点时,常常用到一维搜索,即沿某一已知方向求目标函数 的极小点。一维搜索的方法很多,常用的有:(1)试探法(“成功—失败”,斐波那契法, 0.618 法等);(2)插值法(抛物线插值法,三次插值法等);(3)微积分中的求根法(切 线法,二分法等)。 考虑一维极小化问题 )(min tf bta ≤≤ (2) 若 )(tf 是 ],[ ba 区间上的下单峰函数,我们介绍通过不断地缩短 ],[ ba 的长度,来 搜索得(2)的近似最优解的两个方法。 为了缩短区间 ],[ ba ,逐步搜索得(2)的最优解 *t 的近似值,我们可以采用以下 途径:在 ],[ ba 中任取两个关于 ],[ ba 是对称的点 1t 和 2t (不妨设 12 tt < ,并把它们叫 做搜索点),计算 )( 1tf 和 )( 2tf 并比较它们的大小。对于单峰函数,若 )()( 12 tftf < , 则必有 ],[ 1 * tat ∈ ,因而 ],[ 1ta 是缩短了的单峰区间;若 )()( 21 tftf < ,则有 ],[ 2 * btt ∈ ,故 ],[ 2 bt 是缩短了的单峰区间;若 )()( 12 tftf = ,则 ],[ 1ta 和 ],[ 2 bt 都是 缩短了的单峰。因此通过两个搜索点处目标函数值大小的比较,总可以获得缩短了的单 峰区间。对于新的单峰区间重复上述做法,显然又可获得更短的单峰区间。如此进行, 在单峰区间缩短到充分小时,我们可以取最后的搜索点作为(2)最优解的近似值。 应该按照怎样的规则来选取探索点,使给定的单峰区间的长度能尽快地缩短? 2.1.1 Fibonacci 法 如用 nF 表示计算 n 个函数值能缩短为单位长区间的最大原区间长度,可推出 nF 满 足关系 110 == FF ,,3,2,12 L=+= −− nFFF nnn 数列 }{ nF 称为 Fibonacci 数列, nF 称为第 n 个 Fibonacci 数,相邻两个 Fibonacci 数之 比 n n F F 1− 称为 Fibonacci 分数。 当用斐波那契法以 n 个探索点来缩短某一区间时,区间长度的第一次缩短率为 n n F F 1− ,其后各次分别为 2 1 2 3 1 2 ,,, F F F F F F n n n n L − − − − 。由此,若 1t 和 )( 122 ttt < 是单峰区间 ],[ ba 中第 1 个和第 2 个探索点的话,那么应有比例关系 n n F F ab at 11 −=− − , n n F F ab at 22 −=− − 从而 )(1 1 abF Fat n n −+= − , )(2 2 abF Fat n n −+= − (3) 它们关于 ],[ ba 确是对称的点。 -37- 如果要求经过一系列探索点搜索之后,使最后的探索点和最优解之间的距离不超 过精度 0>δ ,这就要求最后区间的长度不超过δ ,即 δ≤− nF ab (4) 据此,我们应按照预先给定的精度δ ,确定使(4)成立的最小整数 n 作为搜索次数, 直到进行到第 n 个探索点时停止。 用上述不断缩短函数 )(tf 的单峰区间 ],[ ba 的办法,来求得问题(2)的近似解, 是 Kiefer(1953 年)提出的,叫做 Finbonacci 法,具体步骤如下: o1 选取初始数据,确定单峰区间 ],[ 00 ba ,给出搜索精度 0>δ ,由(4)确定搜 索次数 n 。 o2 00 ,,1 bbaak === ,计算最初两个搜索点,按(3)计算 1t 和 2t 。 o3 while 1−< nk )(),( 2211 tfftff == if 21 ff < )()( )1(;; 1122 abknF knFatttta −− −−+=== else )()( )1(;; 2211 baknF knFbttttb −− −−+=== end 1+= kk end o4 当进行至 1−= nk 时, )(2 1 21 batt +== 这就无法借比较函数值 )( 1tf 和 )( 2tf 的大小确定最终区间,为此,取 ⎪⎪ ⎩ ⎪⎪ ⎨ ⎧ −++= += ))(2 1( )(2 1 1 2 abat bat ε 其中ε 为任意小的数。在 1t 和 2t 这两点中,以函数值较小者为近似极小点,相应的函数 值为近似极小值。并得最终区间 ],[ 1ta 或 ],[ 2 bt 。 由上述分析可知,斐波那契法使用对称搜索的方法,逐步缩短所考察的区间,它能 以尽量少的函数求值次数,达到预定的某一缩短率。 例 3 试用斐波那契法求函数 2)( 2 +−= tttf 的近似极小点,要求缩短后的区间 不大于区间 ]3,1[− 的 0.08 倍。 程序留作习题。 2.1.2 0.618 法 若 0>ω ,满足比例关系 -38- ω ωω −= 1 1 称之为黄金分割数,其值为 L6180339887.02 15 =−=ω 。 黄金分割数ω 和 Fibonacci 分数之间有着重要的关系 n n n F F 1lim − ∞→ =ω 。 现用不变的区间缩短率 0.618,代替斐波那契法每次不同的缩短率,就得到了黄金 分割法(0.618 法)。这个方法可以看成是斐波那契法的近似,实现起来比较容易,效果 也相当好,因而易于为人们所接受。 用 0.618 法求解,从第 2 个探索点开始每增加一个探索点作一轮迭代以后,原单 峰区间要缩短 0.618 倍。计算 n 个探索点的函数值可以把原区间 ],[ 00 ba 连续缩短 1−n 次,因为每次的缩短率均为 μ ,故最后的区间长度为 1 00 )( −− nab μ 这就是说,当已知缩短的相对精度为δ 时,可用下式计算探索点个数 n : δμ ≤−1n 当然,也可以不预先计算探索点的数目 n ,而在计算过程中逐次加以判断,看是否已满 足了提出的精度要求。 0.618 法是一种等速对称进行试探的方法,每次的探索点均取在区间长度的 0.618 倍和 0.382 倍处。 2.2 二次插值法 对极小化问题(2),当 )(tf 在 ],[ ba 上连续时,可以考虑用多项式插值来进行一 维搜索。它的基本思想是:在搜索区间中,不断用低次(通常不超过三次)多项式来近 似目标函数,并逐步用插值多项式的极小点来逼近(2)的最优解。 2.3 无约束极值问题的解法 无约束极值问题可表述为 )(),( min nExxf ∈ (5) 求解问题(5)的迭代法大体上分为两点:一是用到函数的一阶导数或二阶导数, 称为解析法。另一是仅用到函数值,称为直接法。 2.3.1 解析法 2.3.1.1 梯度法(最速下降法) 对基本迭代格式 k k kk ptxx +=+1 (6) 我们总是考虑从点 kx 出发沿哪一个方向 kp ,使目标函数 f 下降得最快。微积分的知识 告诉我们,点 kx 的负梯度方向 )( kk xfp −∇= , 是从点 kx 出发使 f 下降最快的方向。为此,称负梯度方向 )( kxf∇− 为 f 在点 kx 处的 最速下降方向。 -39- 按基本迭代格式(6),每一轮从点 kx 出发沿最速下降方向 )( kxf∇− 作一维搜索, 来建立求解无约束极值问题的方法,称之为最速下降法。 这个方法的特点是,每轮的搜索方向都是目标函数在当前点下降最快的方向。同 时,用 0)( =∇ kxf 或 ε≤∇ )( kxf 作为停止条件。其具体步骤如下: 1°选取初始数据。选取初始点 0x ,给定终止误差,令 0:=k 。 2°求梯度向量。计算 )( kxf∇ , 若 ε≤∇ )( kxf ,停止迭代,输出 kx 。否则, 进行 3°。 3° 构造负梯度方向。取 )( kk xfp −∇= . 4° 进行一维搜索。求 kt ,使得 )(min)( 0 kk t k k k tpxfptxf +=+ ≥ 令 ,1:,1 +=+=+ kkptxx k k kk 转 2°。 例 4 用最速下降法求解无约束非线性规划问题 2 2 2 1 25)(min xxxf += 其中 Txxx ),( 21= ,要求选取初始点 Tx )2,2(0 = 。 解 (i) Txxxf )50,2()( 21=∇ 编写 M 文件 detaf.m,定义函数 )(xf 及其梯度列向量如下 function [f,df]=detaf(x); f=x(1)^2+25*x(2)^2; df=[2*x(1) 50*x(2)]; (ii)编写主程序文件zuisu.m如下: 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/2; f=detaf(x+t*p); end x=x+t*p; [f0,g]=detaf(x); end x,f0 2.3.1.2 Newton 法 考虑目标函数 f 在点 kx 处的二次逼近式 ))(()(2 1 )()()()()( 2 kkTkkTkk xxxfxxxxxfxfxQxf −∇−+−∇+=≈ 假定 Hesse 阵 -40- ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∂ ∂ ∂∂ ∂ ∂∂ ∂ ∂ ∂ =∇ 2 2 1 2 1 2 2 1 2 2 )()( )()( )( n k n k n kk k x xf xx xf xx xf x xf xf L MLM L 正定。 由于 )(2 kxf∇ 正定,函数Q 的驻点 1+kx 是 )(xQ 的极小点。为求此极小点,令 0))(()()( 121 =−∇+∇=∇ ++ kkkkk xxxfxfxQ , 即可解得 )()]([ 121 kkkk xfxfxx ∇∇−= −+ . 对照基本迭代格式(1),可知从点 kx 出发沿搜索方向。 )()]([ 12 kkk xfxfp ∇∇−= − 并取步长 1=kt 即可得 )(xQ 的最小点 1+kx 。通常,把方向 kp 叫做从点 kx 出发的 Newton 方向。从一初始点开始,每一轮从当前迭代点出发,沿 Newton 方向并取步长 为 1 的求解方法,称之为 Newton 法。其具体步骤如下: 1°选取初始数据。选取初始点 0x ,给定终止误差 0>ε ,令 0:=k 。 2°求梯度向量。计算 )( kxf∇ ,若 ε≤∇ )( kxf ,停止迭代,输出 kx 。否则, 进行 3°。 3°构造 Newton 方向。计算 12 )]([ −∇ kxf ,取 )()]([ 12 kkk xfxfp ∇∇−= − . 4° 求下一迭代点。令 1:,1 +=+=+ kkpxx kkk ,转 2°。 例 5 用 Newton 法求解, 2 2 2 1 4 2 4 1 25)(min xxxxxf ++= 选取 Tx )2,2(0 = 。 解:(i) Txxxxxxxf ]210024[)( 2 2 1 3 2 2 21 3 1 ++=∇ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + +=∇ 2 1 2 221 21 2 2 2 12 23004 4212 xxxx xxxxf (ii)编写 M 文件 nwfun.m 如下: function [f,df,d2f]=nwfun(x); f=x(1)^4+25*x(2)^4+x(1)^2*x(2)^2; df=[4*x(1)^3+2*x(1)*x(2)^2;100*x(2)^3+2*x(1)^2*x(2)]; d2f=[2*x(1)^2+2*x(2)^2,4*x(1)*x(2) 4*x(1)*x(2),300*x(2)^2+2*x(1)^2]; (III)编写主程序文件 example5.m 如下: clc x=[2;2]; [f0,g1,g2]=nwfun(x); while norm(g1)>0.00001 p=-inv(g2)*g1; x=x+p; -41- [f0,g1,g2]=nwfun(x); end x, f0 如果目标函数是非二次函数,一般地说,用 Newton 法通过有限轮迭代并不能保证 可求得其最优解。 为了提高计算精度,我们在迭代时可以采用变步长计算上述问题,编写主程序文件 example5_2 如下: clc,clear x=[2;2]; [f0,g1,g2]=nwfun(x); while norm(g1)>0.00001 p=-inv(g2)*g1;p=p/norm(p); t=1.0;f=nwfun(x+t*p); while f>f0 t=t/2;f=nwfun(x+t*p); end x=x+t*p; [f0,g1,g2]=nwfun(x); end x,f0 Newton 法的优点是收敛速度快;缺点是有时不好用而需采取改进措施,此外,当 维数较高时,计算 12 )]([ −∇− kxf 的工作量很大。 2.3.1.3 变尺度法 变尺度法(Variable Metric Algorithm)是近 20 多年来发展起来的,它不仅是求解 无约束极值问题非常有效的算法,而且也已被推广用来求解约束极值问题。由于它既避 免了计算二阶导数矩阵及其求逆过程,又比梯度法的收敛速度快,特别是对高维问题具 有显著的优越性,因而使变尺度法获得了很高的声誉。下面我们就来简要地介绍一种变 尺度法—DFP 法的基本原理及其计算过程。这一方法首先由 Davidon 在 1959 年提出, 后经 Fletcher 和 Powell 加以改进。 我们已经知道,牛顿法的搜索方向是 )()]([ 12 kk xfxf ∇∇− − ,为了不计算二阶 导数矩阵 )]([ 2 kxf∇ 及其逆阵,我们设法构造另一个矩阵,用它来逼近二阶导数矩阵 的逆阵 12 )]([ −∇ kxf ,这一类方法也称拟牛顿法(Quasi-Newton Method)。 下面研究如何构造这样的近似矩阵,并将它记为 )(kH 。我们要求:每一步都能以 现有的信息来确定下一个搜索方向;每做一次选代,目标函数值均有所下降;这些近似 矩阵最后应收敛于解点处的 Hesse 阵的逆阵。 当 )(xf 是二次函数时,其 Hesse 阵为常数阵 A ,任两点 kx 和 1+kx 处的梯度之差为 )()()( 11 kkkk xxAxfxf −=∇−∇ ++ 或 )]()([ 111 kkkk xfxfAxx ∇−∇=− +−+ 对于非二次函数,仿照二次函数的情形,要求其 Hesse 阵的逆阵的第 1+k 次近似 矩阵 )1( +kH 满足关系式 )]()([ 1)1(1 kkkkk xfxfHxx ∇−∇=− +++ (7) -42- 这就是常说的拟 Newton 条件。 若令 ⎩ ⎨ ⎧ −=Δ ∇−∇=Δ + + kkk kkk xxx xfxfG 1 1)( )()( (8) 则式(7)变为 )()1( kkk GHx Δ=Δ + , (9) 现假定 )(kH 已知,用下式求 )1( +kH (设 )(kH 和 )1( +kH 均为对称正定阵); )()()1( kkk HHH Δ+=+ (10) 其中 )(kHΔ 称为第 k 次校正矩阵。显然, )1( +kH 应满足拟 Newton 条件(9),即要求 )()()( )( kkkk GHHx ΔΔ+=Δ 或 )()()()( kkkkk GHxGH Δ−Δ=ΔΔ (11) 由此可以设想, )(kHΔ 的一种比较简单的形式是 TkkkTkkk WGHQxH )()( )()()()()( Δ−Δ=Δ (12) 其中 )(kQ 和 )(kW 为两个待定列向量。 将式(12)中的 )(kHΔ 代入(11),得 )()()()()()()()( )()( kkkkTkkkkTkk GHxGWGHGQx Δ−Δ=ΔΔ−ΔΔ 这说明,应使 1)()( )()()()( =Δ=Δ kTkkTk GWGQ (13) 考虑到 )(kHΔ 应为对称阵,最简单的办法就是取 ⎪⎩ ⎪⎨ ⎧ Δ= Δ= )()()( )( kk k k k k k GHW xQ ξ η (14) 由式(13)得 1)()( )()()()( =ΔΔ=ΔΔ kkTk k kTk k GHGGx ξη (15) 若 )()( kTk Gx ΔΔ 和 )()()( )( kkTk GHG ΔΔ 不等于零,则有 ⎪ ⎪ ⎩ ⎪⎪ ⎨ ⎧ ΔΔ= ΔΔ=ΔΔ= )()()( )()( )( 1 )( 1 )( 1 kkTkk kTkkTkk GHG xGGx ξ η (16) 于是,得校正矩阵 )()()( )()()()( )( )( )( )( )( )( kkTk kTkkk kTk Tkk k GHG HGGH xG xxH ΔΔ ΔΔ−ΔΔ ΔΔ=Δ (17) 从而得到 )()()( )()()()( )( )()1( )( )( )( )( kkTk kTkkk kTk Tkk kk GHG HGGH xG xxHH ΔΔ ΔΔ−ΔΔ ΔΔ+=+ (18) 上述矩阵称为尺度矩阵。通常,我们取第一个尺度矩阵 )0(H 为单位阵,以后的尺度矩 -43- 阵按式(18)逐步形成。可以证明: (i)当 kx 不是极小点且 )(kH 正定时,式(17)右端两项的分母不为零,从而可 按式(18)产生下一个尺度矩阵 )1( +kH ; (ii)若 )(kH 为对称正定阵,则由式(18)产生的 )1( +kH 也是对称正定阵; (iii)由此推出 DFP 法的搜索方向为下降方向。 现将 DFP 变尺度法的计算步骤总结如下。 1°给定初始点 0x 及梯度允许误差 0>ε 。 2°若 ε≤∇ )( 0xf ,则 0x 即为近似极小点,停止迭代,否则,转向下一步。 3°令 IH =)0( (单位矩阵), )( 0)0(0 xfHp ∇−= 在 0p 方向进行一维搜索,确定最佳步长 0λ : )()(min 0 0 000 pxfpxf λλ λ +=+ 如此可得下一个近似点 0 0 01 pxx λ+= 4°一般地,设已得到近似点 kx ,算出 )( kxf∇ ,若 ε≤∇ )( 0xf 则 kx 即为所求的近似解,停止迭代;否则,计算 )(kH : )1()1()1( )1()1()1()1( 1)1( 11 )1()( )( )( )( )( −−− −−−− −− −− − ΔΔ ΔΔ−ΔΔ ΔΔ+= kkTk kTkkk kTk Tkk kk GHG HGGH xG xxHH 并令 )()( kkk xfHp ∇−= ,在 kp 方向上进行一维搜索,得 kλ ,从而可得下一个近似点 k k kk pxx λ+=+1 5°若 1+kx 满足精度要求,则 1+kx 即为所求的近似解,否则,转回 4°,直到求 出某点满足精度要求为止。 2.3.2 直接法 在无约束非线性规划方法中,遇到问题的目标函数不可导或导函数的解析式难以 表示时,人们一般需要使用直接搜索方法。同时,由于这些方法一般都比较直观和易于 理解,因而在实际应用中常为人们所采用。下面我们介绍 Powell 方法。 这个方法主要由所谓基本搜索、加速搜索和调整搜索方向三部分组成,具体步骤 如下: 1° 选取初始数据。选取初始点 0x , n 个线性无关初始方向,组成初搜索方向组 },,,{ 110 −nppp L 。给定终止误差 0>ε ,令 0:=k 。 2°进行基本搜索。令 kxy =:0 ,依次沿 },,,{ 110 −nppp L 中的方向进行一维搜 索。对应地得到辅助迭代点 nyyy ,,, 21 L ,即 1 1 1 − − − += j j jj ptyy njtpyfptyf jj t j j j ,,1),(min)( 11 0 1 1 1 L=+=+ −− ≥ − − − -44- 3°构造加速方向。令 0yyp nn −= ,若 ε≤np ,停止迭代,输出 nk yx =+1 。 否则进行 4°。 4°确定调整方向。按下式 }1|)()(max{)()( 11 njyfyfyfyf jjmm ≤≤−=− −− 找出 m 。若 )]()([2)2()(2)( 100 mmnn yfyfyyfyfyf −<−+− − 成立,进行 5°。否则,进行 6°。 5°调整搜索方向组。令 )(min)(: 0 1 nn t n n nn n nk tpyfptyfptyx +=++= ≥ + . 同时,令 },,,,,,{:},,,{ 1110 1 110 nnmm k n pppppppp −+− + + = LLL , 1: += kk ,转 2°。 6°不调整搜索方向组。令 1:,:1 +==+ kkyx nk ,转 2°。 2.4 Matlab 求无约束极值问题 在 Matlab 工具箱中,用于求解无约束极值问题的函数有 fminunc 和 fminsearch,用 法介绍如下。 求函数的极小值 , )(min xf x 其中 x 可以为标量或向量。 Matlab 中 fminunc 的基本命令是 [X,FVAL]=FMINUNC(FUN,X0,OPTIONS,P1,P2, ...) 其中的返回值 X 是所求得的极小点,FVAL 是函数的极小值,其它返回值的含义参见相 关的帮助。FUN 是一个 M 文件,当 FUN 只有一个返回值时,它的返回值是函数 )(xf ; 当 FUN 有两个返回值时,它的第二个返回值是 )(xf 的梯度向量;当 FUN 有三个返回 值时,它的第三个返回值是 )(xf 的二阶导数阵(Hessian 阵)。X0 是向量 x 的初始值, OPTIONS 是优化参数,可以使用缺省参数。P1,P2 是可以传递给 FUN 的一些参数。 例 6 求函数 2 1 22 12 )1()(100)( xxxxf −+−= 的最小值。 解:编写 M 文件 fun2.m 如下: function [f,g]=fun2(x); f=100*(x(2)-x(1)^2)^2+(1-x(1))^2; g=[-400*x(1)*(x(2)-x(1)^2)-2*(1-x(1));200*(x(2)-x(1)^2)]; 编写主函数文件example6.m如下: options = optimset('GradObj','on'); [x,y]=fminunc('fun2',rand(1,2),options) 即可求得函数的极小值。 在求极值时,也可以利用二阶导数,编写 M 文件 fun3.m 如下: function [f,df,d2f]=fun3(x); f=100*(x(2)-x(1)^2)^2+(1-x(1))^2; df=[-400*x(1)*(x(2)-x(1)^2)-2*(1-x(1));200*(x(2)-x(1)^2)]; d2f=[-400*x(2)+1200*x(1)^2+2,-400*x(1) -400*x(1),200]; -45- 编写主函数文件example62.m如下: options = optimset('GradObj','on','Hessian','on'); [x,y]=fminunc('fun3',rand(1,2),options) 即可求得函数的极小值。 求多元函数的极值也可以使用 Matlab 的 fminsearch 命令,其使用格式为: [X,FVAL,EXITFLAG,OUTPUT]=FMINSEARCH(FUN,X0,OPTIONS,P1,P2,...) 例 7 求函数 3)sin()( += xxf 取最小值时的 x 值。 解 编写 )(xf 的 M 文件 fun4.m 如下: function f=fun4(x); f=sin(x)+3; 编写主函数文件example7.m如下: x0=2; [x,y]=fminsearch(@fun4,x0) 即求得在初值 2 附近的极小点及极小值。 §3 约束极值问题 带有约束条件的极值问题称为约束极值问题,也叫规划问题。 求解约束极值问题要比求解无约束极值问题困难得多。为了简化其优化工作,可采 用以下方法:将约束问题化为无约束问题;将非线性规划问题化为线性规划问题,以及 能将复杂问题变换为较简单问题的其它方法。 库恩—塔克条件是非线性规划领域中最重要的理论成果之一,是确定某点为最优点 的必要条件,但一般说它并不是充分条件(对于凸规划,它既是最优点存在的必要条件, 同时也是充分条件)。 3.1 二次规划 若某非线性规划的目标函数为自变量 x 的二次函数,约束条件又全是线性的,就称 这种规划为二次规划。 Matlab 中二次规划的数学模型可表述如下: ,2 1min xfHxx TT + s.t. ⎩ ⎨ ⎧ =⋅ ≤ beqxAeq bAx 这里 H 是实对称矩阵, bf , 是列向量, A 是相应维数的矩阵。 Matlab 中求解二次规划的命令是 [X,FVAL]= QUADPROG(H,f,A,b,Aeq,beq,LB,UB,X0,OPTIONS) 返回值 X 是决策向量 x 的值,返回值 FVAL 是目标函数在 x 处的值。(具体细节可以参 看在 Matlab 指令中运行 help quadprog 后的帮助)。 例 8 求解二次规划 ⎪ ⎪ ⎩ ⎪⎪ ⎨ ⎧ ≥ ≤+ ≤+ −−+−= 0, 94 3 36442)(min 21 21 21 21 2 221 2 1 xx xx xx xxxxxxxf 解 编写如下程序: -46- 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)) 求得 1.02501)( Min, 0500.1 9500.1 −=⎥⎦ ⎤ ⎢⎣ ⎡= xfx 。 3.2 罚函数法 利用罚函数法,可将非线性规划问题的求解,转化为求解一系列无约束极值问题, 因而也称这种方法为序列无约束最小化技术,简记为 SUMT (Sequential Unconstrained Minization Technique)。 罚函数法求解非线性规划问题的思想是,利用问题中的约束函数作出适当的罚函 数,由此构造出带参数的增广目标函数,把问题转化为无约束非线性规划问题。主要有 两种形式,一种叫外罚函数法,另一种叫内罚函数法,下面介绍外罚函数法。 考虑问题: )(min xf s.t. ⎪⎩ ⎪⎨ ⎧ == =≥ =≤ ,,1,0)( s,,1,,0)( ,,,1 ,0)( tmxk jxh rixg m j i L L L 取一个充分大的数 0>M ,构造函数 ∑∑∑ === +−+= t i i s i i r i i xkMxhMxgMxfMxP 111 |)(|)0),(min()0),(max()(),( (或 ||)(||0 )(minsum0 )(maxsum)(),( xKMxHMxGMxfMxP +⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛−⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛+= 这里 [])()()( 1 xgxgxG rL= , [ ])()()( 1 xhxhxH sL= , [])()()( 1 tkxkxK tL= ,Matlab 中可以直接利用 max 、min 和 sum 函数。)则以 增广目标函数 ),( MxP 为目标函数的无约束极值问题 ),(min MxP 的最优解 x 也是原问题的最优解。 例 9 求下列非线性规划 ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ ≥ =+−− ≥− ++= .0, 02 0 8)(min 21 2 21 2 2 1 2 2 2 1 xx xx xx xxxf 解 (i)编写 M 文件 test.m function g=test(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)+... -47- M*abs(-x(1)-x(2)^2+2); 或者是利用Matlab的求矩阵的极小值和极大值函数编写test.m如下: function g=test(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); 我们也可以修改罚函数的定义,编写test.m如下: function g=test(x); M=50000; f=x(1)^2+x(2)^2+8; g=f-M*min(min(x),0)-M*min(x(1)^2-x(2),0)+M*(-x(1)-x(2)^2+2)^2; (ii)在 Matlab 命令窗口输入 [x,y]=fminunc('test',rand(2,1)) 即可求得问题的解。 3.3 Matlab 求约束极值问题 在 Matlab 优化工具箱中,用于求解约束最优化问题的函数有:fminbnd、fmincon、 quadprog、fseminf、fminimax,上面我们已经介绍了函数 fmincon 和 quadprog。 3.3.1 fminbnd 函数 求单变量非线性函数在区间上的极小值 ],[ , )(min 21 xxxxf x ∈ Matlab 的命令为 [X,FVAL] = FMINBND(FUN,x1,x2,OPTIONS), 它的返回值是极小点 x 和函数的极小值。这里 fun 是用 M 文件定义的函数或 Matlab 中 的单变量数学函数。 例 10 求函数 ]5,0[ ,1)3()( 2 ∈−−= xxxf 的最小值。 解 编写 M 文件 fun5.m function f=fun5(x); f=(x-3)^2-1; 在 Matlab 的命令窗口输入 [x,y]=fminbnd('fun5',0,5) 即可求得极小点和极小值。 3.3.2 fseminf 函数 求 }0),(,0)(,0)(|)({min ≤=≤ wxPHIxCeqxCxF x ⎩ ⎨ ⎧ = ≤ BeqxAeq BxA * * s.t. 其中 ),(),(),( wxPHIxCeqxC 都是向量函数; w 是附加的向量变量, w 的每个分量都 限定在某个区间内。 上述问题的 Matlab 命令格式为 X=FSEMINF(FUN,X0,NTHETA,SEMINFCON,A,B,Aeq,Beq) 其中 FUN 用于定义目标函数 )(xF ;X0 为 x 的初始值;NTHETA 是半无穷约束 ),( wxPHI 的个数;函数 SEMINFCON 用于定义非线性不等式约束 )(xC ,非线性等 -48- 式约束 )(xCeq 和半无穷约束 ),( wxPHI 的每一个分量函数,函数 SEMINFCON 有两个 输入参量 X 和 S,S 是推荐的取样步长,也许不被使用。 例 11 求函数 2 3 2 2 2 1 )5.0()5.0()5.0()( −+−+−= xxxxf 取最小值时的 x 值, 约束为: 1)sin()50(1000 1)cos()sin(),( 331 2 1211111 ≤−−−−= xxwwxwxwwxK 1)sin()50(1000 1)cos()sin(),( 332 2 2122222 ≤−−−−= xxwwxwxwwxK 1001 1 ≤≤ w , 1001 2 ≤≤ w 解 (1)编写 M 文件 fun6.m 定义目标函数如下: function f=fun6(x,s); f=sum((x-0.5).^2); (2)编写 M 文件 fun7.m 定义约束条件如下: function [c,ceq,k1,k2,s]=fun7(x,s); c=[];ceq=[]; if isnan(s(1,1)) s=[0.2,0;0.2 0]; end %取样值 w1=1:s(1,1):100; w2=1:s(2,1):100; %半无穷约束 k1=sin(w1*x(1)).*cos(w1*x(2))-1/1000*(w1-50).^2-sin(w1*x(3))-x(3)-1; k2=sin(w2*x(2)).*cos(w2*x(1))-1/1000*(w2-50).^2-sin(w2*x(3))-x(3)-1; %画出半无穷约束的图形 plot(w1,k1,'-',w2,k2,'+'); (3)调用函数 fseminf 在 Matlab 的命令窗口输入 [x,y]=fseminf(@fun6,rand(3,1),2,@fun7) 即可。 3.3.3 fminimax 函数 求解 ⎭⎬⎫ ⎩⎨⎧ )(maxmin xF iFx ⎪ ⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ ≤≤ = ≤ = ≤ UBxLB xCeq xC BeqxAeq bxA 0)( 0)( * * .s.t 其中 )}(,),({)( 1 xFxFxF mL= 。 上述问题的 Matlab 命令为 X=FMINIMAX(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON) -49- 例 12 求函数族 )}(),(),(),(),({ 54321 xfxfxfxfxf 取极大极小值时的 x 值。其中: ⎪ ⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ −+= −−= −+= −−= +−−+= 8)( )( 183)( 3)( 30440482)( 215 214 213 2 2 2 12 21 2 2 2 11 xxxf xxxf xxxf xxxf xxxxxf 解 (1)编写 M 文件 fun8.m 定义向量函数如下: function f=fun8(x); f=[2*x(1)^2+x(2)^2-48*x(1)-40*x(2)+304 -x(1)^2-3*x(2)^2 x(1)+3*x(2)-18 -x(1)-x(2) x(1)+x(2)-8]; (2)调用函数 fminimax [x,y]=fminimax(@fun8,rand(2,1)) 3.3.4 利用梯度求解约束优化问题 例 13 已知函数 )12424()( 221 2 2 2 1 1 ++++= xxxxxexf x ,且满足非线性约束: ⎩ ⎨ ⎧ −≥ −≤−− 10 5.1 21 2121 xx xxxx 求 )(min xf x 分析:当使用梯度求解上述问题时,效率更高并且结果更准确。 题目中目标函数的梯度为: ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ++ +++++ )244( )168424( 21 2121 2 2 2 1 1 1 xxe xxxxxxe x x 解 (1)编写 M 文件 fun9.m 定义目标函数及梯度函数: function [f,df]=fun9(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)]; (2)编写 M 文件 fun10.m 定义约束条件及约束条件的梯度函数: function [c,ceq,dc,dceq]=fun10(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=[]; (3)调用函数 fmincon,编写主函数文件 example13.m 如下: %采用标准算法 options=optimset('largescale','off'); %采用梯度 options=optimset(options,'GradObj','on','GradConstr','on'); [x,y]=fmincon(@fun9,rand(2,1),[],[],[],[],[],[],@fun10,options) -50- 3.4 Matlab 优化工具箱的用户图形界面解法 Matlab 优化工具箱中的 optimtool 命令提供了优化问题的用户图形界面解法。 optimtool 可应用到所有优化问题的求解,计算结果可以输出到 Matlab 工作空间中。 图 1 优化问题用户图形界面解法示意图 例 14 用 optimtool 重新求解例 1。 利用例 1 已经定义好的函数 fun1 和 fun2。在 Matlab 命令窗口运行 optimtool,就 打 开图形界面,如图 1 所示,填入有关的参数,未填入的参数取值为空或者为默认值,然 后用鼠标点一下“start”按钮,就得到求解结果,再使用“file”菜单下的“Export to Workspace…”选项,把计算结果输出到 Matlab 工作空间中去。 §4 飞行管理问题 在约 10,000m 高空的某边长 160km 的正方形区域内,经常有若干架飞机作水平 飞行。区域内每架飞机的位置和速度向量均由计算机记录其数据,以便进行飞行管理。 当一架欲进入该区域的飞机到达区域边缘时,记录其数据后,要立即计算并判断是否会 与区域内的飞机发生碰撞。如果会碰撞,则应计算如何调整各架(包括新进入的)飞机 飞行的方向角,以避免碰撞。现假定条件如下: 1)不碰撞的标准为任意两架飞机的距离大于 8km; 2)飞机飞行方向角调整的幅度不应超过 30 度; 3)所有飞机飞行速度均为每小时 800km; 4)进入该区域的飞机在到达区域边缘时,与区域内飞机的距离应在 60km 以上; 5)最多需考虑 6 架飞机; 6)不必考虑飞机离开此区域后的状况。 请你对这个避免碰撞的飞行管理问题建立数学模型,列出计算步骤,对以下数据进 行计算(方向角误差不超过 0.01 度),要求飞机飞行方向角调整的幅度尽量小。 设该区域 4 个顶点的座标为(0,0),(160,0),(160,160),(0,160)。记录数据见表 1。 -51- 表 1 飞行记录数据 飞机编号 横座标 x 纵座标 y 方向角(度) 1 150 140 243 2 85 85 236 3 150 155 220.5 4 145 50 159 5 130 150 230 新进入 0 0 52 注:方向角指飞行方向与 x 轴正向的夹角。 为方便以后的讨论,我们引进如下记号: D 为飞行管理区域的边长; Ω 为飞行管理区域,取直角坐标系使其为 ],0[],0[ DD × ; a 为飞机飞行速度, 800=a km/h; ),( 00 ii yx 为第i 架飞机的初始位置; ))(),(( tytx ii 为第i 架飞机在t 时刻的位置; 0 iθ 为第i 架飞机的原飞行方向角,即飞行方向与 x 轴夹角, πθ 20 0 <≤ i ; iθΔ 为第i 架飞机的方向角调整, 66 πθπ ≤Δ≤− i ; iii θθθ Δ+= 0 为第i 架飞机调整后的飞行方向角。 4.1 模型一 根据相对运动的观点在考察两架飞机i 和 j 的飞行时,可以将飞机i 视为不动而飞 机 j 以相对速度 )sinsin,coscos( ijijijij aaaavvv θθθθ −−=−= (19) 相对于飞机i 运动,对(19)式进行适当的化约可得 )2cos,2sin(2sin2 ijijijav θθθθθθ ++−−= ))22cos(),22(cos(2sin2 ijijija θθπθθπθθ ++++−= (20) 不妨设 ij θθ ≥ ,此时相对飞行方向角为 22 ji ij θθπβ ++= ,见图 2。 图 2 相对飞行方向角 由于两架飞机的初始距离为 200200 )()()0( jijiij yyxxr −+−= (21) -52- )0( 8arcsin0 ij ij r=α (22) 则只要当相对飞行方向角 ijβ 满足 00 2 ijijij απβα −<< (23) 时,两架飞机不可能碰撞(见图 2)。 记 0 ijβ 为调整前第 j 架飞机相对于第i 架飞机的相对速度(矢量)与这两架飞机连 线(从 j 指向i 的矢量)的夹角(以连线矢量为基准,逆时针方向为正,顺时针方向为 负)。则由式(23)知,两架飞机不碰撞的条件为 00 )(2 1 ijjiij αθθβ >Δ+Δ+ (24) 其中 0 mnβ =相对速度 mnv 的幅角-从 n 指向 m 的连线矢量的幅角 )()(arg nnmm ii iyxiyx ee mn +−+ −= θθ (注意 0 mnβ 表达式中的i 表示虚数单位)这里我们利用复数的幅角,可以很方便地计算 角度 0 mnβ ( 6,,2,1, L=nm )。 本问题中的优化目标函数可以有不同的形式:如使所有飞机的最大调整量最小; 所有飞机的调整量绝对值之和最小等。这里以所有飞机的调整量绝对值之和最小为目标 函数,可以得到如下的数学规划模型: ∑ = Δ 6 1 min i iθ s.t. 00 )(2 1 ijjiij αθθβ >Δ+Δ+ , 6,,2,1, L=ji , ji ≠ o30≤Δ iθ , 6,,2,1 L=i 利用如下的程序: clc,clear x0=[150 85 150 145 130 0]; y0=[140 85 155 50 150 0]; q=[243 236 220.5 159 230 52]; xy0=[x0; y0]; d0=dist(xy0); %求矩阵各个列向量之间的距离 d0(find(d0==0))=inf; a0=asind(8./d0) %以度为单位的反函数 xy1=x0+i*y0 xy2=exp(i*q*pi/180) for m=1:6 for n=1:6 if n~=m b0(m,n)=angle((xy2(n)-xy2(m))/(xy1(m)-xy1(n))); end -53- end end b0=b0*180/pi; dlmwrite('txt1.txt',a0,'delimiter', '\t','newline','PC'); fid=fopen('txt1.txt','a'); fwrite(fid,'~','char'); %往纯文本文件中写 LINGO 数据的分割符 dlmwrite('txt1.txt',b0,'delimiter', '\t','newline','PC','-append','roffset', 1) 求得 0 ijα 的值如表 2 所示。 表 2 0 ijα 的值 1 2 3 4 5 6 1 0 5.39119 32.23095 5.091816 20.96336 2.234507 2 5.39119 0 4.804024 6.61346 5.807866 3.815925 3 32.23095 4.804024 0 4.364672 22.83365 2.125539 4 5.091816 6.61346 4.364672 0 4.537692 2.989819 5 20.96336 5.807866 22.83365 4.537692 0 2.309841 6 2.234507 3.815925 2.125539 2.989819 2.309841 0 求得 0 ijβ 的值如表 3 所示。 表 3 0 ijβ 的值 1 2 3 4 5 6 1 0 109.26 -128.25 24.18 173.07 14.475 2 109.26 0 -88.871 -42.244 -92.305 9 3 -128.25 -88.871 0 12.476 -58.786 0.31081 4 24.18 -42.244 12.476 0 5.9692 -3.5256 5 173.07 -92.305 -58.786 5.9692 0 1.9144 6 14.475 9 0.31081 -3.5256 1.9144 0 上述飞行管理的数学规划模型可如下输入 LINGO 求解: model: sets: plane/1..6/:delta; link(plane,plane):alpha,beta; endsets data: alpha=@file('txt1.txt'); !需要在alpha的数据后面加上分隔符"~"; beta=@file('txt1.txt'); enddata min=@sum(plane:@abs(delta)); @for(plane:@bnd(-30,delta,30)); @for(link(i,j)|i#ne#j:@abs(beta(i,j)+0.5*delta(i)+0.5*delta(j))>a lpha(i,j)); end 求得的最优解为 o83858.23 =Δθ , o0138.215 −=Δθ , o7908.06 =Δθ ,其它调整 角度为 0。 4.2 模型二 -54- 两架飞机 ji, 不发生碰撞的条件为 64))()(())()(( 22 >−+− tytytxtx jiji (25) 51 ≤≤ i , 61 ≤≤+ ji , },min{0 ji TTt ≤≤ 其中 ji TT , 分别表示第 ji, 架飞机飞出正方形区域边界的时刻。这里 iii atxtx θcos)( 0 += , iii atyty θsin)( 0 += , ni ,,2,1 L= ; iii θθθ Δ+= 0 , 6|| πθ ≤Δ i , ni ,,2,1 L= ; 下面我们把约束条件(25)加强为对所有的时间t 都成立,记 ),(~),(~),(~64))()(())()(( 222 , jictjibtjiatytytxtxl jijiji ++=−−+−= 其中 2sin4),(~ 22 jiajia θθ −= , )]sin))(sin0()0(()cos))(cos0()0([(2),(~ jijijiji yyxxajib θθθθ −−+−−= 64))0()0(())0()0((),(~ 22 −−+−= jiji yyxxjic 则两架 ji, 飞机不碰撞的条件是 0),(~),(~4),(~),( 2 <−=Δ jicjiajibji (26) 这样我们建立如下的非线性规划模型 ∑ = Δ 6 1 2)( i iθ s.t. 0),( <Δ ji , 51 ≤≤ i , 61 ≤≤+ ji 6 πθ ≤Δ i , 6,,2,1 L=i 习 题 三 1. 用最速下降法(梯度法)求函数: 2 221 2 121 22264)( xxxxxxxf −−−+= 的极大点。给定初始点 Tx )1,1(0 = 。 2. 试用牛顿法求解: 2 1)(min 2 2 2 1 ++−= xxxf 取初始点 Tx )0,4()0( = ,并将采用变步长和采用固定步长 0.1=λ 时的情形做比较。 3. 某工厂向用户提供发动机,按合同规定,其交货数量和日期是:第一季度末交 40 台,第二季末交 60 台,第三季末交 80 台。工厂的最大生产能力为每季 100 台,每 季的生产费用是 22.050)( xxxf += (元),此处 x 为该季生产发动机的台数。若工厂 生产的多,多余的发动机可移到下季向用户交货,这样,工厂就需支付存贮费,每台发 动机每季的存贮费为 4 元。问该厂每季应生产多少台发动机,才能既满足交货合同,又 使工厂所花费的费用最少(假定第一季度开始时发动机无存货)。 4. 用 Matlab 的非线性规划命令 fmincon 求解飞行管理问题的模型二。 5. 用罚函数法求解飞行管理问题的模型二。 -55- 6. 求下列问题的解 3 2 22 2 11 332)(max xxxxxxf ++++= s.t. 1022 3 2 22 2 11 ≤++++ xxxxx 503 2 22 2 11 ≤−+++ xxxxx 4022 32 2 11 ≤+++ xxxx 23 2 1 =+ xx 12 21 ≥+ xx 不约束321 ,,0 xxx ≥ 7. 图 3 影院剖面示意图 图 3 为影院的剖面示意图,座位的满意程度主要取决于视角α 和仰角 β 。视 角 α 是 观众眼睛到屏幕上、下边缘视线的夹角,α 越大越好;仰角 β 是观众眼睛到屏幕上边缘 视线与水平线的夹角,β 太大使人的头部过分上仰,引起不舒服感,一般要求 β 不超过 o30 。 记影院屏幕高 h ,上边缘距地面高 H ,地板线倾角θ ,第一排和最后一排座位与屏 幕水平距离分别为 d 和 D ,观众平均座高为 c(指眼睛到地面的距离)。已知参数 8.1=h , 5=H , 5.4=d , 19=D , 1.1=c (单位:m)。 (1)地板线倾角 o10=θ ,问最佳座位在什么地方? (2)求地板线倾角(一般不超过 o20 ),使所有观众的平均满意程度最大。 (3)地板线设计成什么形状可以进一步提高观众的满意程度。 -56- 第四章 动态规划 §1 引言 1.1 动态规划的发展及研究内容 动态规划(dynamic programming)是运筹学的一个分支,是求解决策过程(decision process)最优化的数学方法。20 世纪 50 年代初 R. E. Bellman 等人在研究多阶段决策过 程(multistep decision process)的优化问题时,提出了著名的最优性原理(principle of optimality),把多阶段过程转化为一系列单阶段问题,逐个求解,创立了解决这类过程 优化问题的新方法—动态规划。1957 年出版了他的名著《Dynamic Programming》,这 是该领域的第一本著作。 动态规划问世以来,在经济管理、生产调度、工程技术和最优控制等方面得到了广 泛的应用。例如最短路线、库存管理、资源分配、设备更新、排序、装载等问题,用动 态规划方法比用其它方法求解更为方便。 虽然动态规划主要用于求解以时间划分阶段的动态过程的优化问题,但是一些与时 间无关的静态规划(如线性规划、非线性规划),只要人为地引进时间因素,把它视为 多阶段决策过程,也可以用动态规划方法方便地求解。 应指出,动态规划是求解某类问题的一种方法,是考察问题的一种途径,而不是 一种特殊算法(如线性规划是一种算法)。因而,它不象线性规划那样有一个标准的数 学表达式和明确定义的一组规则,而必须对具体问题进行具体分析处理。因此,在学习 时,除了要对基本概念和方法正确理解外,应以丰富的想象力去建立模型,用创造性的 技巧去求解。 例 1 最短路线问题 图 1 是一个线路网,连线上的数字表示两点之间的距离(或费用)。试寻求一条由 A 到G 距离最短(或费用最省)的路线。 图 1 最短路线问题 例 2 生产计划问题 工厂生产某种产品,每单位(千件)的成本为 1(千元),每次开工的固定成本为 3 (千元),工厂每季度的最大生产能力为 6(千件)。经调查,市场对该产品的需求量第 一、二、三、四季度分别为 2,3,2,4(千件)。如果工厂在第一、二季度将全年的需 求都生产出来,自然可以降低成本(少付固定成本费),但是对于第三、四季度才能上 市的产品需付存储费,每季每千件的存储费为 0.5(千元)。还规定年初和年末这种产品 均无库存。试制定一个生产计划,即安排每个季度的产量,使一年的总费用(生产成本 和存储费)最少。 1.2 决策过程的分类 根据过程的时间变量是离散的还是连续的,分为离散时间决策过程(discrete-time -57- decision process)和连续时间决策过程(continuous-time decision process);根据过程的 演变是确定的还是随机的,分为确定性决策过程(deterministic decision process)和随 机性决策过程(stochastic decision process),其中应用最广的是确定性多阶段决策过程。 §2 基本概念、基本方程和计算方法 2.1 动态规划的基本概念和基本方程 一个多阶段决策过程最优化问题的动态规划模型通常包含以下要素。 2.1.1 阶段 阶段(step)是对整个过程的自然划分。通常根据时间顺序或空间顺序特征来划分阶 段,以便按阶段的次序解优化问题。阶段变量一般用 nk ,,2,1 L= 表示。在例 1 中由 A 出发为 1=k ,由 )2,1( =iBi 出发为 2=k ,依此下去从 )2,1( =iFi 出发为 6=k ,共 6=n 个阶段。在例 2 中按照第一、二、三、四季度分为 4,3,2,1=k ,共四个阶段。 2.1.2 状态 状态(state)表示每个阶段开始时过程所处的自然状况。它应能描述过程的特征并 且无后效性,即当某阶段的状态变量给定时,这个阶段以后过程的演变与该阶段以前各 阶段的状态无关。通常还要求状态是直接或间接可以观测的。 描述状态的变量称状态变量(state variable)。变量允许取值的范围称允许状态集合 (set of admissible states)。用 kx 表示第 k 阶段的状态变量,它可以是一个数或一个向量。 用 kX 表示第 k 阶段的允许状态集合。在例 1 中 2x 可取 21,BB ,或将 iB 定义为 )2,1( =ii ,则 12 =x 或 2 ,而 }2,1{2 =X 。 n 个阶段的决策过程有 1+n 个状态变量, 1+nx 表示 nx 演变的结果。在例 1 中 7x 取 G ,或定义为1,即 17 =x 。 根据过程演变的具体情况,状态变量可以是离散的或连续的。为了计算的方便有时 将连续变量离散化;为了分析的方便有时又将离散变量视为连续的。 状态变量简称为状态。 2.1.3 决策 当一个阶段的状态确定后,可以作出各种选择从而演变到下一阶段的某个状态,这 种选择手段称为决策(decision),在最优控制问题中也称为控制(control)。 描述决策的变量称决策变量(decision variable),变量允许取值的范围称允许决策 集合(set of admissible decisions)。用 )( kk xu 表示第 k 阶段处于状态 kx 时的决策变量, 它是 kx 的函数,用 )( kk xU 表示 kx 的允许决策集合。在例 1 中 )( 12 Bu 可取 21,CC 或 3C , 可记作 3,2,1)1(2 =u ,而 }3,2,1{)1(2 =U 。 决策变量简称决策。 2.1.4 策略 决策组成的序列称为策略(policy)。由初始状态 1x 开始的全过程的策略记作 )( 11 xp n ,即 )}(,),(),({)( 221111 nnn xuxuxuxp L= . 由第 k 阶段的状态 kx 开始到终止状态的后部子过程的策略记作 )( kkn xp ,即 )}(,),({)( nnkkkkn xuxuxp L= , 1,,2,1 −= nk L . 类似地,由第 k 到第 j 阶段的子过程的策略记作 -58- )}(,),({)( jjkkkkj xuxuxp L= . 可供选择的策略有一定的范围,称为允许策略集合(set of admissible policies),用 )(),(),( 11 kkjkknn xPxPxP 表示。 2.1.5. 状态转移方程 在确定性过程中,一旦某阶段的状态和决策为已知,下阶段的状态便完全确定。用 状态转移方程(equation of state transition)表示这种演变规律,写作 .,,2,1),,(1 nkuxTx kkkk L==+ (1) 在例 1 中状态转移方程为 )(1 kkk xux =+ 。 2.1.6. 指标函数和最优值函数 指标函数(objective function)是衡量过程优劣的数量指标,它是定义在全过程和所有 后部子过程上的数量函数,用 ),,,,( 11, ++ nkkknk xxuxV L 表示, nk ,,2,1 L= 。指标函 数应具有可分离性,即 nkV , 可表为 nkkk Vux ,1,, + 的函数,记为 )),,,(,,(),,,,( 111,111, ++++++ = nkknkkkknkkknk xuxVuxxxuxV LL ϕ 并且函数 kϕ 对于变量 nkV ,1+ 是严格单调的。 过程在第 j 阶段的阶段指标取决于状态 jx 和决策 ju ,用 ),( jjj uxv 表示。指标函 数由 ),,2,1( njv j L= 组成,常见的形式有: 阶段指标之和,即 ∑ = ++ = n kj jjjnkkknk uxvxxuxV ),(),,,,( 11, L , 阶段指标之积,即 ∏ = ++ = n kj jjjnkkknk uxvxxuxV ),(),,,,( 11, L , 阶段指标之极大(或极小),即 ),((min)max),,,,( 11, jjjnjknkkknk uxvxxuxV ≤≤++ =L . 这些形式下第 k 到第 j 阶段子过程的指标函数为 ),,,( 1, +jkkjk xuxV L 。 根据状态转移方程指标函数 nkV , 还可以表示为状态 kx 和策略 knp 的函数,即 ),(, knknk pxV 。在 kx 给定时指标函数 nkV , 对 knp 的最优值称为最优值函数(optimal value function),记为 )( kk xf ,即 ),(opt)( , )( knknk xPp kk pxVxf kknkn∈ = , 其中 opt 可根据具体情况取 max 或 min 。 2.1.7 最优策略和最优轨线 使指标函数 nkV , 达到最优值的策略是从 k 开始的后部子过程的最优策略,记作 },,{ *** nkkn uup L= 。 * 1np 是全过程的最优策略,简称最优策略(optimal policy)。从初始 状态 )( * 11 xx = 出发,过程按照 * 1np 和状态转移方程演变所经历的状态序列 },,,{ * 1 * 2 * 1 +nxxx L 称最优轨线(optimal trajectory)。 -59- 2.1.8 递归方程 如下方程称为递归方程 ⎪⎩ ⎪⎨ ⎧ =⊗= = ++ ∈ ++ 1,,)},(),({opt)( 10)( 11 )( 11 Lnkxfuxvxf xf kkkkk xUu kk nn kkk 或 (2) 在上述方程中,当⊗ 为加法时取 0)( 11 =++ nn xf ;当 ⊗ 为乘法时,取 1)( 11 =++ nn xf 。动 态规划递归方程是动态规划的最优性原理的基础,即:最优策略的子策略,构成最优子 策略。用状态转移方程(1)和递归方程(2)求解动态规划的过程,是由 1+= nk 逆 推至 1=k ,故这种解法称为逆序解法。当然,对某些动态规划问题,也可采用顺序解 法。这时,状态转移方程和递归方程分别为: nkuxTx kk r kk ,,1),,( 1 L== + , ⎪⎩ ⎪⎨ ⎧ =⊗= = −+ ∈ + ++ nkxfuxvxf xf kkkkk xUu kk k r kk ,,1)},(),({opt)( 10( 11 )( 1 10 11 L 或) 例 3 用 lingo 求解例 1 最短路线问题。 model: Title Dynamic Programming; sets: vertex/A,B1,B2,C1,C2,C3,C4,D1,D2,D3,E1,E2,E3,F1,F2,G/:L; road(vertex,vertex)/A B1,A B2,B1 C1,B1 C2,B1 c3,B2 C2,B2 C3,B2 C4, C1 D1,C1 D2,C2 D1,C2 D2,C3 D2,C3 D3,C4 D2,C4 D3, D1 E1,D1 E2,D2 E2,D2 E3,D3 E2,D3 E3, E1 F1,E1 F2,E2 F1,E2 F2,E3 F1,E3 F2,F1 G,F2 G/:D; endsets data: D=5 3 1 3 6 8 7 6 6 8 3 5 3 3 8 4 2 2 1 2 3 3 3 5 5 2 6 6 4 3; L=0,,,,,,,,,,,,,,,; enddata @for(vertex(i)|i#GT#1:L(i)=@min(road(j,i):L(j)+D(j,i))); end 纵上所述,如果一个问题能用动态规划方法求解,那么,我们可以按下列步骤,首 先建立起动态规划的数学模型: (i)将过程划分成恰当的阶段。 (ii)正确选择状态变量 kx ,使它既能描述过程的状态,又满足无后效性,同时确 定允许状态集合 kX 。 (iii)选择决策变量 ku ,确定允许决策集合 )( kk xU 。 (iv)写出状态转移方程。 (v)确定阶段指标 ),( kkk uxv 及指标函数 knV 的形式(阶段指标之和,阶段指标之 积,阶段指标之极大或极小等)。 (vi)写出基本方程即最优值函数满足的递归方程,以及端点条件。 §3 逆序解法的计算框图 -60- 以自由终端、固定始端、指标函数取和的形式的逆序解法为例给出计算框图,其它 情况容易在这个基础上修改得到。 一般化的自由终端条件为 1,1,11 ,,2,1),()( ++++ == nininn nixxf Lϕ (3) 其中ϕ 为已知。固定始端条件可表示为 }{}{ * 111 xxX == 。 如果状态 kx 和决策 ku 是连续变量,用数值方法求解时需按照精度要求进行离散 化。设状态 kx 的允许集合为 nkninixX kkkik ,,2,1,,,2,1},,,2,1|{ LLL ==== . 决策 )( kiki xu 的允许集合为 nkninjuU kki j kiki ,,2,1,,,2,1},,,2,1|{ )( LLL ==== . 状态转移方程和阶段指标应对 kx 的每个取值 kix 和 kiu 的每个取值 )( j kiu 计算,即 ),( )( j kikikk uxTT = , ),( )( j kikik uxvv = 。最优值函数应对 kx 的每个取值 kix 计算。基本方 程可以表为 .1,2,,,,,2,1,,,2,1 ),(opt)( )),,((),()( )( )( 1 )()( LLL nkninj xfxf uxTfuxvxf kki ki j k j kik j kikikk j kikikki j k === = += + (4) 图 2 解法框图 -61- 按照(3),(4)逆向计算出 )( * 11 xf ,为全过程的最优值。记状态 kix 的最优决策为 )(* kiki xu ,由 * 1x 和 )(* kiki xu 按照状态转移方程计算出最优状态,记作 * kx 。并得到相应的 最优决策,记作 )( ** kk xu 。于是最优策略为 )}(,),(),({ *** 2 * 2 * 1 * 1 nn xuxuxu L 。 算法程序的框图如图 2 所示。 图的左边部分是函数序列的递推计算,可输出全过程最优值 )( * 11 xf ,如果需要还 可以输出后部子过程最优值函数序列 )( kik xf 和最优决策序列 )(* kik xu 。计算过程中存 )( kik xf 是备计算 1−kf 之用,在 1−kf 算完后可用 1−kf 将 kf 替换掉;存 )(* kik xu 是备右边 部分读 )( ** kk xu 之用。 图的右边部分是最优状态和最优决策序列的正向计算,可输出最优策略 )}(,),(),({ *** 2 * 2 * 1 * 1 nn xuxuxu L 和最优轨线 },,,{ ** 2 * 1 nxxx L 。 §4 动态规划与静态规划的关系 动态规划与静态规划(线性和非线性规划等)研究的对象本质上都是在若干约束条 件下的函数极值问题。两种规划在很多情况下原则上可以相互转换。 动态规划可以看作求决策 nuuu ,,, 21 L 使指标函数 ),,,( 2111 nn uuuxV L, 达到最优 (最大或最小)的极值问题,状态转移方程、端点条件以及允许状态集、允许决策集等 是约束条件,原则上可以用非线性规划方法求解。 一些静态规划只要适当引入阶段变量、状态、决策等就可以用动态规划方法求解。 下面用例子说明。 例 4 用动态规划解下列非线性规划 ∑ = n k kk ug 1 )(max ; s.t. ∑ = ≥= n k kk uau 1 0, . 其中 )( kk ug 为任意的已知函数。 解 按变量 ku 的序号划分阶段,看作 n 段决策过程。设状态为 121 ,,, +nxxx L ,取 问题中的变量 nuuu ,,, 21 L 为决策。状态转移方程为 .,,2,1,, 11 nkuxxax kkk L=−== + 取 )( kk ug 为阶段指标,最优值函数的基本方程为(注意到 01 =+nx ) )]()([max)( 110 ++≤≤ += kkkkxukk xfxgxf kk ; 1,2,,1,,0 L−=≤≤ nnkaxk ; 0)0(1 =+nf . 按照逆序解法求出对应于 kx 每个取值的最优决策 )(* kk xu ,计算至 )(1 af 后即可利 用状态转移方程得到最优状态序列 }{ * kx 和最优决策序列 )}({ ** kk xu 。 与静态规划相比,动态规划的优越性在于: (i)能够得到全局最优解。由于约束条件确定的约束集合往往很复杂,即使指标 函数较简单,用非线性规划方法也很难求出全局最优解。而动态规划方法把全过程化为 -62- 一系列结构相似的子问题,每个子问题的变量个数大大减少,约束集合也简单得多,易 于得到全局最优解。特别是对于约束集合、状态转移和指标函数不能用分析形式给出的 优化问题,可以对每个子过程用枚举法求解,而约束条件越多,决策的搜索范围越小, 求解也越容易。对于这类问题,动态规划通常是求全局最优解的唯一方法。 (ii)可以得到一族最优解。与非线性规划只能得到全过程的一个最优解不同,动 态规划得到的是全过程及所有后部子过程的各个状态的一族最优解。有些实际问题需要 这样的解族,即使不需要,它们在分析最优策略和最优值对于状态的稳定性时也是很有 用的。当最优策略由于某些原因不能实现时,这样的解族可以用来寻找次优策略。 (iii)能够利用经验提高求解效率。如果实际问题本身就是动态的,由于动态规划 方法反映了过程逐段演变的前后联系和动态特征,在计算中可以利用实际知识和经验提 高求解效率。如在策略迭代法中,实际经验能够帮助选择较好的初始策略,提高收敛速 度。 动态规划的主要缺点是: (i)没有统一的标准模型,也没有构造模型的通用方法,甚至还没有判断一个问 题能否构造动态规划模型的准则。这样就只能对每类问题进行具体分析,构造具体的模 型。对于较复杂的问题在选择状态、决策、确定状态转移规律等方面需要丰富的想象力 和灵活的技巧性,这就带来了应用上的局限性。 (ii)用数值方法求解时存在维数灾(curse of dimensionality)。若一维状态变量有 m 个取值,那么对于 n 维问题,状态 kx 就有 nm 个值,对于每个状态值都要计算、存储函 数 )( kk xf ,对于 n 稍大的实际问题的计算往往是不现实的。目前还没有克服维数灾的 有效的一般方法。 §5 若干典型问题的动态规划模型 5.1 最短路线问题 对于例 1 一类最短路线问题(shortest Path Problem),阶段按过程的演变划分,状 态由各段的初始位置确定,决策为从各个状态出发的走向,即有 )(1 kkk xux =+ ,阶段 指标为相邻两段状态间的距离 ))(,( kkkk xuxd ,指标函数为阶段指标之和,最优值函数 )( kk xf 是由 kx 出发到终点的最短距离(或最小费用),基本方程为 ;1,,)],())(,([min)( 11)( Lnkxfxuxdxf kkkkkkxukk kk =+= ++ .0)( 11 =++ nn xf 利用这个模型可以算出例 l 的最短路线为 GFEDCAB 22121 ,相应的最短距离为 18。 5.2 生产计划问题 对于例 2一类生产计划问题(Production planning problem),阶段按计划时间自然 划分,状态定义为每阶段开始时的储存量 kx ,决策为每个阶段的产量 ku ,记每个阶段 的需求量(已知量)为 kd ,则状态转移方程为 .,,2,1,0,1 nkxduxx kkkkk L=≥−+=+ (5) 设每阶段开工的固定成本费为 a ,生产单位数量产品的成本费为b ,每阶段单位数量产 品的储存费为 c ,阶段指标为阶段的生产成本和储存费之和,即 ⎩ ⎨ ⎧ >++= 0 0,),( kk kkkk ubuacxuxv (6) -63- 指标函数 knV 为 kv 之和。最优值函数 )( kk xf 为从第 k 段的状态 kx 出发到过程终结的最 小费用,满足 .1,,)],(),([min)( 11 Lnkxfuxvxf kkkkkUukk kk =+= ++∈ 其中允许决策集合 kU 由每阶段的最大生产能力决定。若设过程终结时允许存储量为 0 1+nx ,则终端条件是 .0)( 0 11 =++ nn xf (7) (5)~(7)构成该问题的动态规划模型。 5.3 资源分配问题 一种或几种资源(包括资金)分配给若干用户,或投资于几家企业,以获得最大的 效益。资源分配问题(resource allocating Problem)可以是多阶段决策过程,也可以是 静态规划问题,都能构造动态规划模型求解。下面举例说明。 例 5 机器可以在高、低两种负荷下生产。u 台机器在高负荷下的年产量是 )(ug , 在低负荷下的年产量是 )(uh ,高、低负荷下机器的年损耗率分别是 1a 和 1b ( 10 11 <<< ab )。现有 m 台机器,要安排一个 n 年的负荷分配计划,即每年初决定 多少台机器投入高、低负荷运行,使 n 年的总产量最大。如果进一步假设 uug α=)( , uuh β=)( ( 0>> βα ),即高、低负荷下每台机器的年产量分别为α 和 β ,结果将 有什么特点。 解 年度为阶段变量 nk ,,2,1 L= 。状态 kx 为第 k 年初完好的机器数,决策 ku 为 第 k 年投入高负荷运行的台数。当 kx 或 ku 不是整数时,将小数部分理解为一年中正常 工作时间或投入高负荷运行时间的比例。 机器在高、低负荷下的年完好率分别记为 a 和 b ,则 11 aa −= , 11 bb −= ,有 ba < 。因为第 k 年投入低负荷运行的机器台数为 kk ux − ,所以状态转移方程是 )(1 kkkk uxbaux −+=+ (8) 阶段指标 kv 是第 k 年的产量,有 )()(),( kkkkkk uxhuguxv −+= (9) 指标函数是阶段指标之和,最优值函数 )( kk xf 满足 .1,2,,,0 )],(),([max)( 110 Lnkmx xfuxvxf k kkkkkxukk kk =≤≤ += ++≤≤ (10) 及自由终端条件 .0,0)( 111 mxxf nnn ≤≤= +++ (11) 当 kv 中的 hg, 用较简单的函数表达式给出时,对于每个 k 可以用解析方法求解极 值问题。特别,若 uug α=)( , uuh β=)( ,(10)中 的 )](),([ 1 kkkkk xfuxv ++ 将是 ku 的线性函数,最大值点必在区间 kk xu ≤≤0 的左端点 0=ku 或右端点 kk xu = 取得, 即每年初将完好的机器全部投入低负荷或高负荷运行。 §6 具体的应用实例 例 6 设某工厂有 1000 台机器,生产两种产品 BA、 ,若投入 x 台机器生产 A 产 -64- 品,则纯收入为 x5 ,若投入 y 台机器生产 B 种产品,则纯收入为 y4 ,又知:生产 A 种 产品机器的年折损率为 20%,生产 B 产品机器的年折损率为 10%,问在 5 年内如何安 排各年度的生产计划,才能使总收入最高? 解 年度为阶段变量 5,4,3,2,1=k 。 令 kx 表示第 k 年初完好机器数, ku 表示第 k 年安排生产 A 种产品的机器数,则 kk ux − 为第 k 年安排生产 B 种产品的机器数,且 kk xu ≤≤0 。 则第 1+k 年初完好的机器数 kkkkkk uxuxux 1.09.0))(1.01()2.01(1 −=−−+−=+ (12) 令 ),( kkk uxv 表示第 k 年的纯收入, )( kk xf 表示第 k 年初往后各年的最大利润之 和。 显然 0)( 66 =xf (13) 则 )}(),({max)( 110 ++≤≤ += kkkkkxukk xfuxvxf kk )}(4{max)}()(45{max 110110 ++≤≤++≤≤ ++=+−+= kkkkxukkkkkxu xfxuxfuxu kkkk (14) (1) 5=k 时,由(13)、(14)式得 }4{max)( 55055 55 xuxf xu += ≤≤ 55 4xu + 关于 5u 求导,知其导数大于零,所以 55 4xu + 在 5u 等于 5x 处取得最大值, 即 55 xu = 时, 555 5)( xxf = 。 (2) 4=k 时,由(12)、(14)式得 }54{max)( 544044 44 xxuxf xu ++= ≤≤ }5.85.0{max)}1.09.0(54{max 44044440 4444 xuuxxu xuxu +=−++= ≤≤≤≤ 当 44 xu = 时, 444 9)( xxf = (3) 3=k 时, }94{max)( 433033 33 xxuxf xu ++= ≤≤ }1.121.0{max)}1.09.0(94{max 33033330 3333 xuuxxu xuxu +=−++= ≤≤≤≤ 当 33 xu = 时, 333 2.12)( xxf = (4) 2=k 时, }98.1422.0{max}2.124{max)( 220322022 2222 xuxxuxf xuxu +−=++= ≤≤≤≤ 当 02 =u 时, 222 98.14)( xxf = 。 (5) 1=k 时, }482.17498.0{max}98.144{max)( 110211011 1111 xuxxuxf xuxu +−=++= ≤≤≤≤ 当 01 =u 时, 111 482.17)( xxf = 。因为 10001 =x (台) -65- 所以由(12)式,进行回代得 9001.09.0 112 =−= uxx (台) 8101.09.0 223 =−= uxx (台) 6481.09.0 334 =−= uxx (台) 4.5181.09.0 445 =−= uxx (台) 注: 4.5185 =x 台中的 0.4 台应理解为有一台机器只能使用 0.4 年将报废。 例 7 求解下面问题 3 2 21max uuuz = ⎩ ⎨ ⎧ =≥ >=++ 3,2,10 )0(321 iu ccuuu i 解: 按问题的变量个数划分阶段,把它看作为一个三阶段决策问题。设状态变量 为 4321 ,,, xxxx ,并 记 cx =1 ;取问题中的变量 321 ,, uuu 为决策变量;各阶段指标函数 按乘积方式结合。令最优值函数 )( kk xf 表示第 k 阶段的初始状态为 kx ,从 k 阶段到 3 阶段所得到的最大值。 设 33 ux = , 223 xux =+ , cxux ==+ 112 则有 33 xu = , 220 xu ≤≤ , 110 xu ≤≤ 用逆推解法,从后向前依次有 3333 }{max)( 33 xuxf xu == = 及最优解 3 * 3 xu = ),(max)}({max)}({max)( 222022 2 2033 2 2022 222222 xuhuxuxfuxf xuxuxu ≤≤≤≤≤≤ =−== 由 032 2 222 2 2 =−= uxudu dh ,得 22 3 2 xu = 和 02 =u (舍去) 又 222 2 2 2 62 uxdu hd −= ,而 02 2 3 22 2 2 2 22 <−= = x du hd xu ,故 22 3 2 xu = 为极大值点。 所以 3 222 27 4)( xxf = 及最优解 2 * 2 3 2 xu = 。 })(27 4{max)}({max)( 3 1110221011 1111 uxuxfuxf xuxu −== ≤≤≤≤ 同样利用微分法易知 4 111 64 1)( xxf = ,最优解 1 * 1 4 1 xu = 。 由于 1x 已知,因而按计算的顺序反推算,可得各阶段的最优决策和最优值。即 cu 4 1* 1 = , 4 11 64 1)( cxf = 由 cccuxx 4 3 4 1* 112 =−=−= -66- 所以 cxu 2 1 3 2 2 * 2 == , 3 22 16 1)( cxf = 由 cccuxx 4 1 2 1 4 3* 223 =−=−= 所以 cu 4 1* 3 = , cxf 4 1)( 33 = 因此得到最优解为: cu 4 1* 1 = , cu 2 1* 2 = , cu 4 1* 3 = ; 最大值为: 4 1 64 1)(max ccfz == 。 习 题 四 1. 用 Matlab 编程求例 6 的解。 2. 有四个工人,要指派他们分别完成 4 项工作,每人做各项工作所消耗的时间如 表 1 所示。 表 1 工作 工人 A B C D 甲 15 18 21 24 乙 19 23 22 18 丙 26 17 16 19 丁 19 21 23 17 问指派哪个人去完成哪项工作,可使总的消耗时间为最小?试对此问题用动态规划 方法求解。 3. 为保证某一设备的正常运转,需备有三种不同的零件 321 ,, EEE 。若增加备用零 件的数量,可提高设备正常运转的可靠性,但增加了费用,而投资额仅为 8000 元。已 知备用零件数与它的可靠性和费用的关系如表 2 所示。 表 2 增加的可靠性 设备的费用(千元) 备件数 1E 2E 3E 1E 2E 3E 1 2 3 0.3 0.4 0.5 0.2 0.5 0.9 0.1 0.2 0.7 1 2 3 3 5 6 2 3 4 现要求在既不超出投资额的限制,又能尽量提高设备运转的可靠性的条件下,问 各种零件的备件数量应是多少为好? 4. 某工厂购进 100 台机器,准备生产 I、II 两种产品,若生产产品 I,每台机器每 年可收入 45 万元,损坏率为 65%;若生产产品 II,每台机器每年收入为 35 万元,损 坏率为 35%,估计三年后将有新型机器出现,旧的机器将全部淘汰。试问每年应如何 -67- 安排生产,使在三年内收入最多? 5.3 名商人各带 1 名随从乘船渡河,一只小船只能容纳 2 人,由他们自己划行。 随从们密约,在河的任一岸,一旦随从人数比商人多,就杀商人。此密约被商人知道, 如何乘船渡河的大权掌握在商人们手中,商人们怎样安排每次乘船方案,才能安全渡河 呢? 6.某一印刷厂有六项加工任务,对印刷车间和装订车间所需时间(单位:天)如 表 3 所示,试求最优的加工顺序和总加工天数。 表 3 任务 车间 1 2 3 4 5 6 印刷车间 装订车间 3 8 10 12 5 9 2 6 9 5 11 2 -68- 第五章 图与网络模型及方法 §1 概论 图论起源于 18 世纪。第一篇图论论文是瑞士数学家欧拉于 1736 年发表的“哥尼 斯堡的七座桥”。1847 年,克希霍夫为了给出电网络方程而引进了“树”的概念。1857 年,凯莱在计数烷 22 +nn HC 的同分异构物时,也发现了“树”。哈密尔顿于 1859 年提 出“周游世界”游戏,用图论的术语,就是如何找出一个连通图中的生成圈、近几十年 来,由于计算机技术和科学的飞速发展,大大地促进了图论研究和应用,图论的理论和 方法已经渗透到物理、化学、通讯科学、建筑学、运筹学,生物遗传学、心理学、经济 学、社会学等学科中。 图论中所谓的“图”是指某类具体事物和这些事物之间的联系。如果我们用点表示 这些具体事物,用连接两点的线段(直的或曲的)表示两个事物的特定的联系,就得到 了描述这个“图”的几何形象。图论为任何一个包含了一种二元关系的离散系统提供了 一个数学模型,借助于图论的概念、理论和方法,可以对该模型求解。哥尼斯堡七桥问 题就是一个典型的例子。在哥尼斯堡有七座桥将普莱格尔河中的两个岛及岛与河岸联结 起来,问题是要从这四块陆地中的任何一块开始通过每一座桥正好一次,再回到起点。 图 1 哥尼斯堡七桥问题 当然可以通过试验去尝试解决这个问题,但该城居民的任何尝试均未成功。欧拉为了解 决这个问题,采用了建立数学模型的方法。他将每一块陆地用一个点来代替,将每一座 桥用连接相应两点的一条线来代替,从而得到一个有四个“点”,七条“线”的“图”。 问题成为从任一点出发一笔画出七条线再回到起点。欧拉考察了一般一笔画的结构特 点,给出了一笔画的一个判定法则:这个图是连通的,且每个点都与偶数线相关联,将 这个判定法则应用于七桥问题,得到了“不可能走通”的结果,不但彻底解决了这个问 题,而且开创了图论研究的先河。 图与网络是运筹学(Operations Research)中的一个经典和重要的分支,所研究的 问题涉及经济管理、工业工程、交通运输、计算机科学与信息技术、通讯与网络技术等 诸多领域。下面将要讨论的最短路问题、最大流问题、最小费用流问题和匹配问题等都 是图与网络的基本问题。 我们首先通过一些例子来了解网络优化问题。 例 1 最短路问题(SPP-shortest path problem) 一名货柜车司机奉命在最短的时间内将一车货物从甲地运往乙地。从甲地到乙地的 公路网纵横交错,因此有多种行车路线,这名司机应选择哪条线路呢?假设货柜车的运 行速度是恒定的,那么这一问题相当于需要找到一条从甲地到乙地的最短路。 例 2 公路连接问题 某一地区有若干个主要城市,现准备修建高速公路把这些城市连接起来,使得从其 -69- 中任何一个城市都可以经高速公路直接或间接到达另一个城市。假定已经知道了任意两 个城市之间修建高速公路的成本,那么应如何决定在哪些城市间修建高速公路,使得总 成本最小? 例 3 指派问题(assignment problem) 一家公司经理准备安排 N 名员工去完成 N 项任务,每人一项。由于各员工的特点 不同,不同的员工去完成同一项任务时所获得的回报是不同的。如何分配工作方案可以 使总回报最大? 例 4 中国邮递员问题(CPP-chinese postman problem) 一名邮递员负责投递某个街区的邮件。如何为他(她)设计一条最短的投递路线(从 邮局出发,经过投递区内每条街道至少一次,最后返回邮局)?由于这一问题是我国管 梅谷教授 1960 年首先提出的,所以国际上称之为中国邮递员问题。 例 5 旅行商问题(TSP-traveling salesman problem) 一名推销员准备前往若干城市推销产品。如何为他(她)设计一条最短的旅行路线 (从驻地出发,经过每个城市恰好一次,最后返回驻地)?这一问题的研究历史十分悠 久,通常称之为旅行商问题。 例 6 运输问题(transportation problem) 某种原材料有 M 个产地,现在需要将原材料从产地运往 N 个使用这些原材料的工 厂。假定 M 个产地的产量和 N 家工厂的需要量已知,单位产品从任一产地到任一工厂 的运费已知,那么如何安排运输方案可以使总运输成本最低? 上述问题有两个共同的特点:一是它们的目的都是从若干可能的安排或方案中寻求 某种意义下的最优安排或方案,数学上把这种问题称为最优化或优化(optimization) 问题;二是它们都易于用图形的形式直观地描述和表达,数学上把这种与图相关的结构 称为网络(network)。与图和网络相关的最优化问题就是网络最优化或称网络优化 (netwok optimization)问题。所以上面例子中介绍的问题都是网络优化问题。由于多 数网络优化问题是以网络上的流(flow)为研究的对象,因此网络优化又常常被称为网 络流(network flows)或网络流规划等。 下面首先简要介绍图与网络的一些基本概念。 §2 图与网络的基本概念 2.1 无向图 一个无向图(undirected graph)G 是由一个非空有限集合 )(GV 和 )(GV 中某些元素 的无序对集合 )(GE 构成的二元组,记为 ))(),(( GEGVG = 。其中 },,,{)( 21 nvvvGV L= 称为图G 的顶点集(vertex set)或节点集(node set), )(GV 中 的每一个元素 ),,2,1( nivi L= 称为该图的一个顶点(vertex)或节点(node); },,,{)( 21 meeeGE L= 称为图G 的边集(edge set), )(GE 中的每一个元素 ke (即 )(GV 中某两个元素 ji vv , 的无序对) 记为 ),( jik vve = 或 ijjik vvvve == ),,2,1( mk L= , 被称为该图的一条从 iv 到 jv 的边(edge)。 当边 jik vve = 时,称 ji vv , 为边 ke 的端点,并称 jv 与 iv 相邻(adjacent);边 ke 称 为与顶点 ji vv , 关联(incident)。如果某两条边至少有一个公共端点,则称这两条边在 图G 中相邻。 边上赋权的无向图称为赋权无向图或无向网络(undirected network)。我们对图和 网络不作严格区分,因为任何图总是可以赋权的。 -70- 一个图称为有限图,如果它的顶点集和边集都有限。图G 的顶点数用符号 ||V 或 )(Gν 表示,边数用 || E 或 )(Gε 表示。 当讨论的图只有一个时,总是用G 来表示这个图。从而在图论符号中我们常略去 字母G ,例如,分别用 ν,,EV 和ε 代替 )(),(),( GGEGV ν 和 )(Gε 。 端点重合为一点的边称为环(loop)。 一个图称为简单图(simple graph),如果它既没有环也没有两条边连接同一对顶点。 2.2 有向图 定义 一个有向图(directed graph 或 digraph)G 是由一个非空有限集合V 和V 中 某些元素的有序对集合 A 构成的二元组,记为 ),( AVG = 。其中 },,,{ 21 nvvvV L= 称 为图G 的顶点集或节点集, V 中的每一个元素 ),,2,1( nivi L= 称为该图的一个顶点 或节点; },,,{ 21 maaaA L= 称为图G 的弧集(arc set),A 中的每一个元素 ka (即V 中 某两个元素 ji vv , 的有序对) 记为 ),( jik vva = 或 ),,2,1( nkvva jik L== ,被称为该图 的一条从 iv 到 jv 的弧(arc)。 当弧 jik vva = 时,称 iv 为 ka 的尾(tail), jv 为 ka 的头(head),并称弧 ka 为 iv 的 出弧(outgoing arc),为 jv 的入弧(incoming arc)。 对应于每个有向图 D ,可以在相同顶点集上作一个图G ,使得对于 D 的每条弧, G 有一条有相同端点的边与之相对应。这个图称为 D 的基础图。反之,给定任意图G , 对于它的每个边,给其端点指定一个顺序,从而确定一条弧,由此得到一个有向图,这 样的有向图称为G 的一个定向图。 以下若未指明“有向图”三字,“图”字皆指无向图。 2.3 完全图、二分图 每一对不同的顶点都有一条边相连的简单图称为完全图(complete graph)。n 个顶点 的完全图记为 nK 。 若 YXGV U=)( , Φ=YX I , 0|||| ≠YX (这里 || X 表示集合 X 中的元素个 数), X 中无相邻顶点对,Y 中亦然,则称G 为二分图(bipartite graph);特别地,若 YyXx ∈∀∈∀ , ,则 )(GExy ∈ ,则称G 为完全二分图,记成 |||,| YXK 。 2.4 子图 图 H 叫做图G 的 子图(subgraph ),记作 GH ⊂ ,如果 )()( GVHV ⊂ , )()( GEHE ⊂ 。若 H 是G 的子图,则G 称为 H 的母图。 G 的支撑子图(spanning subgraph,又成生成子图)是指满足 )()( GVHV = 的子 图 H 。 2.5 顶点的度 设 )(GVv ∈ ,G 中与 v 关联的边数(每个环算作两条边)称为 v 的度(degree),记 作 )(vd 。若 )(vd 是奇数,称 v 是奇顶点(odd point); )(vd 是偶数,称 v 是偶顶点(even point)。关于顶点的度,我们有如下结果: (i) ∑ ∈ = Vv vd ε2)( (ii) 任意一个图的奇顶点的个数是偶数。 2.6 图与网络的数据结构 -71- 网络优化研究的是网络上的各种优化模型与算法。为了在计算机上实现网络优化的 算法,首先我们必须有一种方法(即数据结构)在计算机上来描述图与网络。一般来说, 算法的好坏与网络的具体表示方法,以及中间结果的操作方案是有关系的。这里我们介 绍计算机上用来描述图与网络的 5 种常用表示方法:邻接矩阵表示法、关联矩阵表示法、 弧表表示法、邻接表表示法和星形表示法。在下面数据结构的讨论中,我们首先假设 ),( AVG = 是一个简单有向图, mAnV == ||,|| ,并假设V 中的顶点用自然数 n,,2,1 L 表示或编号, A 中的弧用自然数 m,,2,1 L 表示或编号。对于有多重边或无向网络的情 况,我们只是在讨论完简单有向图的表示方法之后,给出一些说明。 (i)邻接矩阵表示法 邻接矩阵表示法是将图以邻接矩阵(adjacency matrix)的形式存储在计算机中。图 ),( AVG = 的邻接矩阵是如下定义的:C 是一个 nn × 的 10 − 矩阵,即 nn nnijcC × × ∈= }1,0{)( , ⎩ ⎨ ⎧ ∉ ∈= .),(,0 ,),(,1 Aji Ajicij 也就是说,如果两节点之间有一条弧,则邻接矩阵中对应的元素为 1;否则为 0。 可以看出,这种表示法非常简单、直接。但是,在邻接矩阵的所有 2n 个元素中,只有 m 个为非零元。如果网络比较稀疏,这种表示法浪费大量的存储空间,从而增加了在网络 中查找弧的时间。 图 2 有向图 例 7 对于图 2 所示的有向图,可以用邻接矩阵表示为 ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ 01100 10100 00010 01000 00110 同样,对于网络中的权,也可以用类似邻接矩阵的 nn × 矩阵表示。只是此时一条 弧所对应的元素不再是 1,而是相应的权而已。如果网络中每条弧赋有多种权,则可以 用多个矩阵表示这些权。 (ii)关联矩阵表示法 关联矩阵表示法是将图以关联矩阵(incidence matrix)的形式存储在计算机中.图 ),( AVG = 的关联矩阵 B 是如下定义的: B 是一个 mn × 的矩阵,即 mn mnikbB × × −∈= }1,0,1{)( , -72- ⎪⎩ ⎪⎨ ⎧ ∈=∈∃− ∈=∈∃ = .,0 ,),(, ,1 ,),(,,1 其它 AijkVj AjikVj bik 也就是说,在关联矩阵中,每行对应于图的一个节点,每列对应于图的一条弧。如 果一个节点是一条弧的起点,则关联矩阵中对应的元素为 1;如果一个节点是一条弧的 终点,则关联矩阵中对应的元素为 1− ;如果一个节点与一条弧不关联,则关联矩阵中 对应的元素为 0。对于简单图,关联矩阵每列只含有两个非零元(一个 1+ ,一个 1− )。 可以看出,这种表示法也非常简单、直接。但是,在关联矩阵的所有 nm 个元素中,只 有 m2 个为非零元。如果网络比较稀疏,这种表示法也会浪费大量的存储空间。但由于 关联矩阵有许多特别重要的理论性质,因此它在网络优化中是非常重要的概念。 例 8 对于例 7 所示的图,如果关联矩阵中每列对应弧的顺序为(1,2),(1,3),(2,4), (3,2),(4,3),(4,5),(5,3)和(5,4),则关联矩阵表示为 ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − −− −−− −− 11100000 10110100 01011010 00001101 00000011 同样,对于网络中的权,也可以通过对关联矩阵的扩展来表示。例如,如果网络中 每条弧有一个权,我们可以把关联矩阵增加一行,把每一条弧所对应的权存储在增加的 行中。如果网络中每条弧赋有多个权,我们可以把关联矩阵增加相应的行数,把每一条 弧所对应的权存储在增加的行中。 (iii)弧表表示法 弧表表示法将图以弧表(arc list)的形式存储在计算机中。所谓图的弧表,也就是 图的弧集合中的所有有序对。弧表表示法直接列出所有弧的起点和终点,共需 m2 个存 储单元,因此当网络比较稀疏时比较方便。此外,对于网络图中每条弧上的权,也要对 应地用额外的存储单元表示。例如,例 7 所示的图,假设弧(1,2),(1,3),(2,4),(3,2), (4,3),(4,5),(5,3)和(5,4)上的权分别为 8,9,6,4,0,3,6 和 7,则弧表表示如表 1 所示。 表 1 起点 1 1 2 3 4 4 5 5 终点 2 3 4 2 3 5 3 4 权 8 9 6 4 0 3 6 7 为了便于检索,一般按照起点、终点的字典序顺序存储弧表,如上面的弧表就是按 照这样的顺序存储的。 (iv)邻接表表示法 邻接表表示法将图以邻接表(adjacency lists)的形式存储在计算机中。所谓图的 邻接表,也就是图的所有节点的邻接表的集合;而对每个节点,它的邻接表就是它的所 有出弧。邻接表表示法就是对图的每个节点,用一个单向链表列出从该节点出发的所有 弧,链表中每个单元对应于一条出弧。为了记录弧上的权,链表中每个单元除列出弧的 另一个端点外,还可以包含弧上的权等作为数据域。图的整个邻接表可以用一个指针数 组表示。例如,例 7 所示的图,邻接表表示为 -73- 这是一个 5 维指针数组,每一维(上面表示法中的每一行)对应于一个节点的邻接 表,如第 1 行对应于第 1 个节点的邻接表(即第 1 个节点的所有出弧)。每个指针单元 的第 1 个数据域表示弧的另一个端点(弧的头),后面的数据域表示对应弧上的权。如 第 1 行中的“2”表示弧的另一个端点为 2(即弧为(1,2)),“8”表示对应弧(1,2)上的 权为 8;“3”表示弧的另一个端点为 3(即弧为(1,3)),“9”表示对应弧(1,3)上的权 为 9。又如,第 5 行说明节点 5 出发的弧有(5,3)、(5,4),他们对应的权分别为 6 和 7。 对于有向图 ),( AVG = ,一般用 )(iA 表示节点i 的邻接表,即节点i 的所有出弧构 成的集合或链表(实际上只需要列出弧的另一个端点,即弧的头)。例如上面例子, }3,2{)1( =A , }4,3{)5( =A 等。 (v)星形表示法 星形(star)表示法的思想与邻接表表示法的思想有一定的相似之处。对每个节点, 它也是记录从该节点出发的所有弧,但它不是采用单向链表而是采用一个单一的数组表 示。也就是说,在该数组中首先存放从节点 1 出发的所有弧,然后接着存放从节点 2 出发的所有孤,依此类推,最后存放从节点 n 出发的所有孤。对每条弧,要依次存放其 起点、终点、权的数值等有关信息。这实际上相当于对所有弧给出了一个顺序和编号, 只是从同一节点出发的弧的顺序可以任意排列。此外,为了能够快速检索从每个节点出 发的所有弧,我们一般还用一个数组记录每个节点出发的弧的起始地址(即弧的编号)。 在这种表示法中,可以快速检索从每个节点出发的所有弧,这种星形表示法称为前向星 形(forward star)表示法。 例如,在例 7 所示的图中,仍然假设弧(1,2),(l,3),(2,4),(3,2),(4,3),(4,5), (5,3)和(5,4)上的权分别为 8,9,6,4,0,3,6 和 7。此时该网络图可以用前向 星形表示法表示为表 2 和表 3 。 表 2 节点对应的出弧的起始地址编号数组 节点号 i 1 2 3 4 5 6 起始地址 )(ipoint 1 3 4 5 7 9 表 3 记录弧信息的数组 弧编号 1 2 3 4 5 6 7 8 起点 1 1 2 3 4 4 5 5 终点 2 3 4 2 3 5 3 4 权 8 9 6 4 0 3 6 7 在数组 point 中,其元素个数比图的节点数多 1(即 1+n ),且一定有 1)1( =point , 1)1( +=+ mnpoint 。对于节点i ,其对应的出弧存放在弧信息数组的位置区间为 ]1)1(),([ −+ipointipoint , 如果 )1()( += ipointipoint ,则节点i 没有出弧。这种表示法与弧表表示法也非常相 -74- 似,“记录弧信息的数组”实际上相当于有序存放的“弧表”。只是在前向星形表示法中, 弧被编号后有序存放,并增加一个数组( point )记录每个节点出发的弧的起始编号。 前向星形表示法有利于快速检索每个节点的所有出弧,但不能快速检索每个节点的 所有入弧。为了能够快速检索每个节点的所有入孤,可以采用反向星形(reverse star) 表示法:首先存放进入节点 1 的所有孤,然后接着存放进入节点 2 的所有弧,依此类推, 最后存放进入节点 n 的所有孤。对每条弧,仍然依次存放其起点、终点、权的数值等有 关信息。同样,为了能够快速检索从每个节点的所有入弧,我们一般还用一个数组记录 每个节点的入孤的起始地址(即弧的编号)。例如,例 7 所示的图,可以用反向星形表 示法表示为表 4 和表 5。 表 4 节点对应的入弧的起始地址编号数组 节点号 i 1 2 3 4 5 6 起始地址 )(irpoint 1 1 3 6 8 9 表 5 记录弧信息的数组 弧编号 1 2 3 4 5 6 7 8 终点 2 2 3 3 3 4 4 5 起点 3 1 1 4 5 5 2 4 权 4 8 9 0 6 7 6 3 如果既希望快速检索每个节点的所有出弧,也希望快速检索每个节点的所有入弧, 则可以综合采用前向和反向星形表示法。当然,将孤信息存放两次是没有必要的,可以 只用一个数组(trace)记录一条弧在两种表示法中的对应关系即可。例如,可以在采用 前向星形表示法的基础上,加上上面介绍的 rpoint 数组和如下的trace 数组即可。这相 当于一种紧凑的双向星形表示法,如表 6 所示。 表 6 两种表示法中的弧的对应关系 反向法中弧编号 j 1 2 3 4 5 6 7 8 正向法中弧编号 )( jtrace 4 1 2 5 7 8 3 6 对于网络图的表示法,我们作如下说明: ① 星形表示法和邻接表表示法在实际算法实现中都是经常采用的。星形表示法的 优点是占用的存储空间较少,并且对那些不提供指针类型的语言(如 FORTRAN 语言 等)也容易实现。邻接表表示法对那些提供指针类型的语言(如 C 语言等)是方便的, 且增加或删除一条弧所需的计算工作量很少,而这一操作在星形表示法中所需的计算工 作量较大(需要花费 )(mO 的计算时间)。有关“计算时间”的观念是网络优化中需要 考虑的一个关键因素。 ② 当网络不是简单图,而是具有平行弧(即多重弧)时,显然此时邻接矩阵表示 法是不能采用的。其他方法则可以很方便地推广到可以处理平行弧的情形。 ③ 上述方法可以很方便地推广到可以处理无向图的情形,但由于无向图中边没有 方向,因此可能需要做一些自然的修改。例如,可以在计算机中只存储邻接矩阵的一半 信息(如上三角部分),因为此时邻接矩阵是对称矩阵。无向图的关联矩阵只含有元素 0 和 1+ ,而不含有 1− ,因为此时不区分边的起点和终点。又如,在邻接表和星形表示 法中,每条边会被存储两次,而且反向星形表示显然是没有必要的,等等。 2.7 轨与连通 kk veevevW L2110= ,其 中 )(GEei ∈ , ki ≤≤1 , )(GVv j ∈ , kj ≤≤0 , ie 与 -75- ii vv ,1− 关联,称W 是图G 的一条道路(walk),k 为路长,顶点 0v 和 kv 分别称为W 的起 点和终点,而 121 ,,, −kvvv L 称为它的内部顶点。 若道路W 的边互不相同,则W 称为迹(trail)。若道路W 的顶点互不相同,则W 称 为轨(path)。 称一条道路是闭的,如果它有正的长且起点和终点相同。起点和终点重合的轨叫做 圈(cycle)。 若图G 的两个顶点 vu, 间存在道路,则称u 和 v 连通(connected)。 vu, 间的最短轨 的长叫做 vu, 间的距离。记作 ),( vud 。若图G 的任二顶点均连通,则称G 是连通图。 显然有: (i) 图 P 是一条轨的充要条件是 P 是连通的,且有两个一度的顶点,其余顶点的度 为 2; (ii) 图C 是一个圈的充要条件是C 是各顶点的度均为 2 的连通图。 §3 应用—最短路问题 3.1 两个指定顶点之间的最短路径 问题如下:给出了一个连接若干个城镇的铁路网络,在这个网络的两个指定城镇间, 找一条最短铁路线。 以各城镇为图G 的顶点,两城镇间的直通铁路为图G 相应两顶点间的边,得图G 。 对G 的每一边 e ,赋以一个实数 )(ew —直通铁路的长度,称为 e 的权,得到赋权图G 。 G 的子图的权是指子图的各边的权和。问题就是求赋权图G 中指定的两个顶点 00 ,vu 间的具最小权的轨。这条轨叫做 00 ,vu 间的最短路,它的权叫做 00 ,vu 间的距离,亦记 作 ),( 00 vud 。 求最短路已有成熟的算法:迪克斯特拉(Dijkstra)算法,其基本思想是按距 0u 从 近到远为顺序,依次求得 0u 到G 的各顶点的最短路和距离,直至 0v (或直至G 的所有 顶点),算法结束。为避免重复并保留每一步的计算信息,采用了标号算法。下面是该 算法。 (i) 令 0)( 0 =ul ,对 0uv ≠ ,令 ∞=)(vl , }{ 00 uS = , 0=i 。 (ii) 对每个 iSv ∈ ( ii SVS \= ),用 )}()(),({min uvwulvl iSu + ∈ 代替 )(vl 。计算 )}({min vl iSv∈ ,把达到这个最小值的一个顶点记为 1+iu ,令 }{ 11 ++ = iii uSS U 。 (iii). 若 1|| −= Vi ,停止;若 1|| −< Vi ,用 1+i 代替i ,转(ii)。 算法结束时,从 0u 到各顶点 v 的距离由 v 的最后一次的标号 )(vl 给出。在 v 进入 iS 之前的标号 )(vl 叫 T 标号, v 进入 iS 时的标号 )(vl 叫 P 标号。算法就是不断修改各顶 点的 T 标号,直至获得 P 标号。若在算法运行过程中,将每一顶点获得 P 标号所由来 的边在图上标明,则算法结束时, 0u 至各项点的最短路也在图上标示出来了。 例 9 某公司在六个城市 621 ,,, ccc L 中有分公司,从 ic 到 jc 的直接航程票价记在 下述矩阵的 ),( ji 位置上。( ∞ 表示无直接航路),请帮助该公司设计一张城市 1c 到其它 -76- 城市间的票价最便宜的路线图。 ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∞ ∞ ∞∞ ∞ ∞ 055252510 550102025 25100102040 2010015 252015050 102540500 用矩阵 nna × ( n 为顶点个数)存放各边权的邻接矩阵,行向量 pb 、 1index 、 2index 、 d 分别用来存放 P 标号信息、标号顶点顺序、标号顶点索引、最短通路的值。其中分 量 ⎩ ⎨ ⎧= 顶点未标号当第 顶点已标号当第 i iipb 0 1)( ; )(2 iindex 存放始点到第i 点最短通路中第i 顶点前一顶点的序号; )(id 存放由始点到第i 点最短通路的值。 求第一个城市到其它城市的最短路径的 Matlab 程序如下: clc,clear a=zeros(6); a(1,2)=50;a(1,4)=40;a(1,5)=25;a(1,6)=10; a(2,3)=15;a(2,4)=20;a(2,6)=25; a(3,4)=10;a(3,5)=20; a(4,5)=10;a(4,6)=25; a(5,6)=55; a=a+a'; a(find(a==0))=inf; pb(1:length(a))=0;pb(1)=1;index1=1;index2=ones(1,length(a)); d(1:length(a))=inf;d(1)=0;temp=1; while sum(pb)a(i,k)+a(k,j) a(i,j)=a(i,k)+a(k,j); path(i,j)=k; end end end end a, path 我们使用LINGO9.0编写的FLOYD算法如下: model: sets: -80- nodes/c1..c6/; link(nodes,nodes):w,path; !path标志最短路径上走过的顶点; endsets data: path=0; w=0; @text(mydata1.txt)=@writefor(nodes(i):@writefor(nodes(j): @format(w(i,j),' 10.0f')),@newline(1)); @text(mydata1.txt)=@write(@newline(1)); @text(mydata1.txt)=@writefor(nodes(i):@writefor(nodes(j): @format(path(i,j),' 10.0f')),@newline(1)); enddata calc: w(1,2)=50;w(1,4)=40;w(1,5)=25;w(1,6)=10; w(2,3)=15;w(2,4)=20;w(2,6)=25; w(3,4)=10;w(3,5)=20; w(4,5)=10;w(4,6)=25;w(5,6)=55; @for(link(i,j):w(i,j)=w(i,j)+w(j,i)); @for(link(i,j) |i#ne#j:w(i,j)=@if(w(i,j)#eq#0,10000,w(i,j))); @for(nodes(k):@for(nodes(i):@for(nodes(j): tm=@smin(w(i,j),w(i,k)+w(k,j)); path(i,j)=@if(w(i,j)#gt# tm,k,path(i,j));w(i,j)=tm))); endcalc end §4 树 4.1 基本概念 连通的无圈图叫做树,记之为T 。若图G 满足 )()( TVGV = , )()( GETE ⊂ , 则称T 是G 的生成树。图G 连通的充分必要条件为G 有生成树。一个连通图的生成树 的个数很多,用 )(Gτ 表示G 的生成树的个数,则有公式 公式 (Caylay) 2)( −= n n nKτ 。 公式 )()()( eGeGG ⋅+−= τττ 。 其中 eG − 表示从G 上删除边 e , eG ⋅ 表示把 e 的长度收缩为零得到的图。 树有下面常用的五个充要条件。 定理 1 (i)G 是树当且仅当G 中任二顶点之间有且仅有一条轨道。 (ii)G 是树当且仅当G 无圈,且 1−=νε 。 (iii)G 是树当且仅当G 连通,且 1−=νε 。 (iv)G 是树当且仅当G 连通,且 )(GEe∈∀ , eG − 不连通。 (v)G 是树当且仅当G 无圈, )(GEe∉∀ , eG + 恰有一个圈。 4.2 应用—连线问题 欲修筑连接 n 个城市的铁路,已知i 城与 j 城之间的铁路造价为 ijC ,设计一个线 路图,使总造价最低。 连线问题的数学模型是在连通赋权图上求权最小的生成树。赋权图的具最小权的生 成树叫做最小生成树。 下面介绍构造最小生成树的两种常用算法。 4.2.1 prim 算法构造最小生成树 设置两个集合 P 和Q ,其中 P 用于存放G 的最小生成树中的顶点,集合Q 存放G -81- 的最小生成树中的边。令集合 P 的初值为 }{ 1vP = (假设构造最小生成树时,从顶点 1v 出发),集合Q 的初值为 Φ=Q 。prim 算法的思想是,从所有 Pp ∈ , PVv −∈ 的边 中,选取具有最小权值的边 pv ,将顶点 v 加入集合 P 中,将边 pv 加入集合Q 中,如 此不断重复,直到 VP = 时,最小生成树构造完毕,这时集合Q 中包含了最小生成树 的所有边。 prim 算法如下: (i) }{ 1vP = , Φ=Q ; (ii)while VP =~ 找最小边 pv ,其中 PVvPp −∈∈ , }{vPP += }{pvQQ += end 图 5 最小生成树问题 例 13 用 prim 算法求图 5 的最小生成树。 我们用 nresult ×3 的第一、二、三行分别表示生成树边的起点、终点、权集合。Matlab 程序如下: clc;clear; a=zeros(7); a(1,2)=50; a(1,3)=60; a(2,4)=65; a(2,5)=40; a(3,4)=52;a(3,7)=45; a(4,5)=50; a(4,6)=30;a(4,7)=42; a(5,6)=70; a=a+a';a(find(a==0))=inf; result=[];p=1;tb=2:length(a); while length(result)~=length(a)-1 temp=a(p,tb);temp=temp(:); d=min(temp); [jb,kb]=find(a(p,tb)==d); j=p(jb(1));k=tb(kb(1)); result=[result,[j;k;d]];p=[p,k];tb(find(tb==k))=[]; end result 4.2.1 Kruskal 算法构造最小生成树 科茹斯克尔(Kruskal)算法是一个好算法。Kruskal 算法如下: (i)选 )(1 GEe ∈ ,使得 min)( 1 =ew 。 -82- (ii)若 ieee ,,, 21 L 已选好,则从 },,,{)( 21 ieeeGE L− 中选取 1+ie ,使得 ① }],,,,[{ 121 +ii eeeeG L 中无圈,且 ② min)( 1 =+iew 。 (iii)直到选得 1−νe 为止。 例 14 用 Kruskal 算法构造例 3 的最小生成树。 我们用 nindex ×2 存放各边端点的信息,当选中某一边之后,就将此边对应的顶点序 号中较大序号u 改记为此边的另一序号 v ,同时把后面边中所有序号为u 的改记为 v 。 此方法的几何意义是:将序号u 的这个顶点收缩到 v 顶点,u 顶点不复存在。后面继续 寻查时,发现某边的两个顶点序号相同时,认为已被收缩掉,失去了被选取的资格。 Matlab 程序如下: clc;clear; a(1,2)=50; a(1,3)=60; a(2,4)=65; a(2,5)=40; a(3,4)=52;a(3,7)=45; a(4,5)=50; a(4,6)=30; a(4,7)=42; a(5,6)=70; [i,j,b]=find(a); data=[i';j';b'];index=data(1:2,:); loop=max(size(a))-1; result=[]; while length(result) 的对集 'M ,则 M 称为最大对集;若G 中有一轨,其边交替地在对集 M 内外出现,则 称此轨为 M 的交错轨,交错轨的起止顶点都未被许配时,此交错轨称为可增广轨。 若把可增广轨上在 M 外的边纳入对集,把 M 内的边从对集中删除,则被许配的 顶点数增加 2,对集中的“对儿”增加一个。 1957 年,贝尔热(Berge)得到最大对集的充要条件: 定理 2 M 是图G 中的最大对集当且仅当G 中无 M 可增广轨。 1935 年,霍尔(Hall)得到下面的许配定理: 定理 3 G 为二分图, X 与Y 是顶点集的划分,G 中存在把 X 中顶点皆许配的 -83- 对集的充要条件是, XS ⊂∀ ,则 |||)(| SSN ≥ ,其中 )(SN 是 S 中顶点的邻集。 由上述定理可以得出: 推论 1:若G 是 k 次( )0>k 正则 2 分图,则G 有完美对集。 所谓 k 次正则图,即每顶点皆 k 度的图。 由此推论得出下面的婚配定理: 定理 4 每个姑娘都结识 )1( ≥kk 位小伙子,每个小伙子都结识 k 位姑娘,则每位 姑娘都能和她认识的一个小伙子结婚,并且每位小伙子也能和他认识的一个姑娘结婚。 人员分派问题等实际问题可以化成对集来解决。 人员分派问题:工作人员 nxxx ,,, 21 L 去做 n 件工作 nyyy ,,, 21 L ,每人适合做其 中一件或几件,问能否每人都有一份适合的工作?如果不能,最多几人可以有适合的工 作? 这个问题的数学模型是:G 是二分图,顶点集划分为 YXGV U=)( , },,{ 1 nxxX L= , },,{ 1 nyyY L= ,当且仅当 ix 适合做工作 jy 时, )(GEyx ji ∈ ,求 G 中的最大对集。 解决这个问题可以利用 1965 年埃德门兹(Edmonds)提出的匈牙利算法。 匈牙利算法: (i)从G 中任意取定一个初始对集 M 。 (ii)若 M 把 X 中的顶点皆许配,停止,M 即完美对集;否则取 X 中未被 M 许 配的一顶点u ,记 }{uS = , Φ=T 。 (iii)若 TSN =)( ,停止,无完美对集;否则取 TSNy −∈ )( 。 (iv)若 y 是被 M 许配的,设 Myz ∈ , }{zSS U= , }{yTT U= ,转(iii); 否则,取可增广轨 ),( yuP ,令 ))(())(( MPEPEMM −−= U ,转(ii)。 把以上算法稍加修改就能够用来求二分图的最大完美对集。 最优分派问题:在人员分派问题中,工作人员适合做的各项工作当中,效益未必一 致,我们需要制定一个分派方案,使公司总效益最大。 这个问题的数学模型是:在人员分派问题的模型中,图G 的每边加了权 0)( ≥ji yxw ,表示 ix 干 jy 工作的效益,求加权图G 上的权最大的完美对集。 解决这个问题可以用库恩—曼克莱斯(Kuhn-Munkres)算法。为此,我们要引入 可行顶点标号与相等子图的概念。 定义 若映射 RGVl →)(: ,满足 YyXx ∈∈∀ , , ),()()( yxwylxl ≥+ , 则称l 是二分图G 的可行顶点标号。令 )}()()(),(|{ xywylxlGExyxyEl =+∈= , 称以 lE 为边集的G 的生成子图为相等子图,记作 lG 。 可行顶点标号是存在的。例如 ;),(max)( Xxxywxl Yy ∈= ∈ Yyyl ∈= ,0)( 。 定理 5 lG 的完美对集即为G 的权最大的完美对集。 Kuhn-Munkres 算法 -84- (i)选定初始可行顶点标号l ,确定 lG ,在 lG 中选取一个对集 M 。 (ii)若 X 中顶点皆被 M 许配,停止,M 即G 的权最大的完美对集;否则,取 lG 中未被 M 许配的顶点u ,令 }{uS = , Φ=T 。 (iii)若 TSN lG ⊃)( ,转(iv);若 TSN lG =)( ,取 )}()()({min , xywylxl TySxl −+= ∉∈ α , ⎪⎩ ⎪⎨ ⎧ ∈+ ∈− = 其它),( ,)( ,)( )( vl Tvvl Svvl vl l l α α , ll = , ll GG = 。 (iv)选 TSN lG −)( 中一顶点 y ,若 y 已被 M 许配,且 Myz ∈ ,则 }{zSS U= , }{yTT U= ,转(iii);否则,取 lG 中一个 M 可增广轨 ),( yuP ,令 ))(())(( MPEPEMM −−= U , 转(ii)。 其中 )(SN lG 是 lG 中 S 的相邻顶点集。 §6 Euler 图和 Hamilton 图 6.1 基本概念 定义 经过G 的每条边的迹叫做G 的 Euler 迹;闭的 Euler 迹叫做 Euler 回路或 E 回路;含 Euler 回路的图叫做 Euler 图。 直观地讲,Euler 图就是从一顶点出发每边恰通过一次能回到出发点的那种图,即 不重复地行遍所有的边再回到出发点。 定理 6 (i)G 是 Euler 图的充分必要条件是G 连通且每顶点皆偶次。 ( ii ) G 是 Euler 图的充分必要条件是G 连通且 U d i iCG 1= = , iC 是圈, )()()( jiCECE ji ≠Φ=I 。 (iii)G 中有 Euler 迹的充要条件是G 连通且至多有两个奇次点。 定义 包含G 的每个顶点的轨叫做 Hamilton(哈密顿)轨;闭的 Hamilton 轨叫做 Hamilton 圈或 H 圈;含 Hamilton 圈的图叫做 Hamilton 图。 直观地讲,Hamilton 图就是从一顶点出发每顶点恰通过一次能回到出发点的那种 图,即不重复地行遍所有的顶点再回到出发点。 6.2 Euler 回路的 Fleury 算法 1921 年,Fleury 给出下面的求 Euler 回路的算法。 Fleury 算法: 1o. )(0 GVv ∈∀ ,令 00 vW = 。 2o. 假设迹 iii vevevW L110= 已经选定,那么按下述方法从 },,{ 1 ieeE L− 中选取 边 1+ie : (i) 1+ie 和 iv 相关联; -85- (ii)除非没有别的边可选择,否则 1+ie 不是 },,{ 1 ii eeGG L−= 的割边(cut edge)。 (所谓割边是一条删除后使连通图不再连通的边)。 3o. 当第 2 步不能再执行时,算法停止。 6.3 应用 6.3.1 邮递员问题 中国邮递员问题 一位邮递员从邮局选好邮件去投递,然后返回邮局,当然他必须经过他负责投递的 每条街道至少一次,为他设计一条投递路线,使得他行程最短。 上述中国邮递员问题的数学模型是:在一个赋权连通图上求一个含所有边的回路, 且使此回路的权最小。 显然,若此连通赋权图是 Euler 图,则可用 Fleury 算法求 Euler 回路,此回路即为 所求。 对于非 Euler 图,1973 年,Edmonds 和 Johnson 给出下面的解法: 设G 是连通赋权图。 (i)求 )}2(mod1)(),(|{0 =∈= vdGVvvV 。 (ii)对每对顶点 0, Vvu ∈ ,求 ),( vud ( ),( vud 是u 与 v 的距离,可用 Floyd 算法 求得)。 (iii)构造完全赋权图 || 0VK ,以 0V 为顶点集,以 ),( vud 为边uv 的权。 (iv)求 || 0VK 中权之和最小的完美对集 M 。 (v)求 M 中边的端点之间的在G 中的最短轨。 (vi)在(v)中求得的每条最短轨上每条边添加一条等权的所谓“倍边”(即共端 点共权的边)。 (vii)在(vi)中得的图 'G 上求 Euler 回路即为中国邮递员问题的解。 多邮递员问题: 邮局有 )2( ≥kk 位投递员,同时投递信件,全城街道都要投递,完成任务返回邮 局,如何分配投递路线,使得完成投递任务的时间最早?我们把这一问题记成 kPP。 kPP 的数学模型如下: ),( EVG 是连通图, )(0 GVv ∈ ,求G 的回路 kCC ,,1 L ,使得 (i) )(0 iCVv ∈ , ki ,,2,1 L= , (ii) min)(max )(1 =∑ ∈≤≤ iCEeki ew , (iii)U k i i GECE 1 )()( = = 6.3.2 旅行商(TSP)问题 一名推销员准备前往若干城市推销产品,然后回到他的出发地。如何为他设计一条 最短的旅行路线(从驻地出发,经过每个城市恰好一次,最后返回驻地)?这个问题称 为旅行商问题。用图论的术语说,就是在一个赋权完全图中,找出一个有最小权的 Hamilton 圈。称这种圈为最优圈。与最短路问题及连线问题相反,目前还没有求解旅行 商问题的有效算法。所以希望有一个方法以获得相当好(但不一定最优)的解。 一个可行的办法是首先求一个 Hamilton 圈C ,然后适当修改C 以得到具有较小权 的另一个 Hamilton 圈。修改的方法叫做改良圈算法。设初始圈 121 vvvvC nL= 。 -86- (i)对于 nji <<+< 11 ,构造新的 Hamilton 圈: 12112121 vvvvvvvvvvvC njjijjjiij LLL +++−−= , 它是由C 中删去边 1+iivv 和 1+jjvv ,添加边 jivv 和 11 ++ ji vv 而得到的。若 )()()()( 1111 ++++ +<+ jjiijiji vvwvvwvvwvvw ,则 以 ijC 代替C , ijC 叫做C 的改良圈。 (ii)转(i),直至无法改进,停止。 用改良圈算法得到的结果几乎可以肯定不是最优的。为了得到更高的精确度,可以 选择不同的初始圈,重复进行几次算法,以求得较精确的结果。 这个算法的优劣程度有时能用 Kruskal 算法加以说明。假设C 是G 中的最优圈。 则对于任何顶点 v , vC − 是在 vG − 中的 Hamilton 轨,因而也是 vG − 的生成树。由 此推知:若T 是 vG − 中的最优树,同时 e 和 f 是和 v 关联的两条边,并使得 )()( fwew + 尽可能小,则 )()()( fwewTw ++ 将是 )(Cw 的一个上界。 这里介绍的方法已被进一步发展。圈的修改过程一次替换三条边比一次仅替换两条 边更为有效;然而,有点奇怪的是,进一步推广这一想法,就不对了。 例 15 从北京(Pe)乘飞机到东京(T)、纽约(N)、墨西哥城(M)、伦敦(L)、巴黎(Pa) 五城市做旅游,每城市恰去一次再回北京,应如何安排旅游线,使旅程最短?各城市之 间的航线距离如表 7。 表 7 六城市间的距离 L M N Pa Pe T L 56 35 21 51 60 M 56 21 57 78 70 N 35 21 36 68 68 Pa 21 57 36 51 61 Pe 51 78 68 51 13 T 60 70 68 61 13 解:编写程序如下: function main clc,clear global a a=zeros(6); a(1,2)=56;a(1,3)=35;a(1,4)=21;a(1,5)=51;a(1,6)=60; a(2,3)=21;a(2,4)=57;a(2,5)=78;a(2,6)=70; a(3,4)=36;a(3,5)=68;a(3,6)=68; a(4,5)=51;a(4,6)=61; a(5,6)=13; a=a+a'; L=size(a,1); c1=[5 1:4 6]; [circle,long]=modifycircle(c1,L); c2=[5 6 1:4];%改变初始圈,该算法的最后一个顶点不动 [circle2,long2]=modifycircle(c2,L); if long20 flag=0; for m=1:L-3 for n=m+2:L-1 if a(c1(m),c1(n))+a(c1(m+1),c1(n+1))<... a(c1(m),c1(m+1))+a(c1(n),c1(n+1)) flag=1; c1(m+1:n)=c1(n:-1:m+1); end end end end long=a(c1(1),c1(L)); for i=1:L-1 long=long+a(c1(i),c1(i+1)); end circle=c1; 6.3.3 旅行商问题的数学表达式 设城市的个数为 n , ijd 是两个城市i 与 j 之间的距离, 0=ijx 或 1(1 表示走过城 市i 到城市 j 的路,0 表示没有选择走这条路)。则有 ∑ ≠ ji ijij xdmin s.t. 1 1 =∑ = n j ijx , ni ,,2,1 L= ,(每个点只有一条边出去) 1 1 =∑ = n i ijx , nj ,,2,1 L= ,(每个点只有一条边进去) 1|| , −≤∑ ∈ sx sji ij , 1||2 −≤≤ ns , },,2,1{ ns L⊂ (除起点和终点外,各边不构成圈) }1,0{∈ijx , nji ,,2,1, L= , ji ≠ 。 其中 || s 表示集合 s 中元素个数。 将旅行商问题写成数学规划的具体形式还需要一定的技巧,下面的例子我们引用 LINGO 帮助中的一个程序。 例 16 已知 SV 地区各城镇之间距离见表 8,某公司计划在 SV 地区做广告宣传, 推销员从城市 1 出发,经过各个城镇,再回到城市 1。为节约开支,公司希望推销员走 过这 10 个城镇的总距离最少。 表 8 城镇之间的距离 2 3 4 5 6 7 8 9 10 1 8 5 9 12 14 12 16 17 22 2 9 15 17 8 11 18 14 22 3 7 9 11 7 12 12 17 4 3 17 10 7 15 18 -88- 5 8 10 6 15 15 6 9 14 8 16 7 8 6 11 8 11 11 9 10 解 编写 LINGO 程序如下: MODEL: SETS: CITY / 1.. 10/: U; ! U( I) = sequence no. of city; LINK( CITY, CITY): DIST, ! The distance matrix; X; ! X( I, J) = 1 if we use link I, J; ENDSETS DATA: !Distance matrix, it need not be symmetric; DIST =0 8 5 9 12 14 12 16 17 22 8 0 9 15 17 8 11 18 14 22 5 9 0 7 9 11 7 12 12 17 9 15 7 0 3 17 10 7 15 18 12 17 9 3 0 8 10 6 15 15 14 8 11 17 8 0 9 14 8 16 12 11 7 10 10 9 0 8 6 11 16 18 12 7 6 14 8 0 11 11 17 14 12 15 15 8 6 11 0 10 22 22 17 18 15 16 11 11 10 0; ENDDATA !The model:Ref. Desrochers & Laporte, OR Letters, Feb. 91; N = @SIZE( CITY); MIN = @SUM( LINK: DIST * X); @FOR( CITY( K): ! It must be entered; @SUM( CITY( I)| I #NE# K: X( I, K)) = 1; ! It must be departed; @SUM( CITY( J)| J #NE# K: X( K, J)) = 1; ! Weak form of the subtour breaking constraints; ! These are not very powerful for large problems; @FOR( CITY( J)| J #GT# 1 #AND# J #NE# K: U( J) >= U( K) + X ( K, J) - ( N - 2) * ( 1 - X( K, J)) + ( N - 3) * X( J, K))); ! Make the X's 0/1; @FOR( LINK: @BIN( X)); ! For the first and last stop we know...; @FOR( CITY( K)| K #GT# 1: U( K) <= N - 1 - ( N - 2) * X( 1, K); U( K) >= 1 + ( N - 2) * X( K, 1)); END §7 最大流问题 7.1 最大流问题的数学描述 7.1.1 网络中的流 定义 在以V 为节点集, A 为弧集的有向图 ),( AVG = 上定义如下的权函数: -89- (i) RAL →: 为孤上的权函数,弧 Aji ∈),( 对应的权 ),( jiL 记为 ijl ,称为孤 ),( ji 的容量下界(lower bound); (ii) RAU →: 为弧上的权函数,弧 Aji ∈),( 对应的权 ),( jiU 记为 iju ,称为孤 ),( ji 的容量上界,或直接称为容量(capacity); (iii) RVD →: 为顶点上的权函数,节点 Vi ∈ 对应的权 )(iD 记为 id ,称为顶 点i 的供需量(supply/demand); 此时所构成的网络称为流网络,可以记为 ),,,,( DULAVN = 。 由于我们只讨论 AV, 为有限集合的情况,所以对于弧上的权函数 UL, 和顶点上的 权函数 D ,可以直接用所有孤上对应的权和顶点上的权组成的有限维向量表示,因此 DUL ,, 有时直接称为权向量,或简称权。由于给定有向图 ),( AVG = 后,我们总是可 以在它的弧集合和顶点集合上定义各种权函数,所以流网络一般也直接简称为网络。 在流网络中,弧 ),( ji 的容量下界 ijl 和容量上界 iju 表示的物理意义分别是:通过该 弧发送某种“物质”时,必须发送的最小数量为 ijl ,而发送的最大数量为 iju 。顶点 Vi ∈ 对应的供需量 id 则表示该顶点从网络外部获得的“物质”数量( 0>id 时),或从该顶 点发送到网络外部的“物质”数量( 0id 时,表示有 id 个单位的流量从网络外部流入该顶点,因此顶点i 称 为供应点(supply node)或源(source),有时也形象地称为起始点或发点等;当 0ijf ,则称这条轨 P 为从 s 到t 的关于 f 的可增广轨。 令 ⎩ ⎨ ⎧ − = 为后向弧当 为前向弧当 ( ( ),, ),, jif jifu ij ijij ijδ , }min{ ijδδ = 则在这条可增广轨上每条前向弧的流都可以增加一个量δ ,而相应的后向弧的流可减 少δ ,这样就可使得网络的流量获得增加,同时可以使每条弧的流量不超过它的容量, 而且保持为正,也不影响其它弧的流量。总之,网络中 f 可增广轨的存在是有意义的, 因为这意味着 f 不是最大流。 7.3 最大流的一种算法—标号法 标号法是由 Ford 和 Fulkerson 在 1957 年提出的。用标号法寻求网络中最大流的基 本思想是寻找可增广轨,使网络的流量得到增加,直到最大为止。即首先给出一个初始 流,这样的流是存在的,例如零流。如果存在关于它的可增广轨,那么调整该轨上每条 弧上的流量,就可以得到新的流。对于新的流,如果仍存在可增广轨,则用同样的方法 使流的值增大,继续这个过程,直到网络中不存在关于新得到流的可增广轨为止,则该 流就是所求的最大流。 这种方法分为以下两个过程: A.标号过程:通过标号过程寻找一条可增广轨。 B.增流过程:沿着可增广轨增加网络的流量。 这两个过程的步骤分述如下。 -92- (A)标号过程: (i)给发点标号为 ),( ∞+s 。 (ii)若顶点 x 已经标号,则对 x 的所有未标号的邻接顶点 y 按以下规则标号: ① 若 Ayx ∈),( ,且 xyxy uf < 时,令 },min{ xxyxyy fu δδ −= , 则给顶点 y 标号为 ),( yx δ+ ,若 xyxy uf = ,则不给顶点 y 标号。 ② Axy ∈),( ,且 0>yxf ,令 },min{ xyxy f δδ = ,则给 y 标号为 ),( yx δ− ,若 0=yxf ,则不给 y 标号。 (iii)不断地重复步骤(ii)直到收点t 被标号,或不再有顶点可以标号为止。当t 被标号时,表明存在一条从 s 到t 的可增广轨,则转向增流过程(B)。如若t 点不能被 标号,且不存在其它可以标号的顶点时,表明不存在从 s 到t 的可增广轨,算法结束, 此时所获得的流就是最大流。 (B)增流过程 (i)令 tu = 。 (ii)若u 的标号为 tv δ,( + ),则 tvuvu ff δ+= ;若u 的标号为 ),( tv δ− ,则 tuvuv ff δ−= 。 (iii)若 su = ,把全部标号去掉,并回到标号过程(A)。否则,令 vu = ,并回 到增流过程(ii)。 求网络 ),,,,( UAVtsN = 中的最大流 x 的算法的程序设计具体步骤如下: 对每个节点 j ,其标号包括两部分信息 f(j))max),(pred( j 该节点在可能的增广路中的前一个节点 )(pred j ,以及沿该可能的增广路到该节点为止 可以增广的最大流量 )(fmax j 。 STEP0 置初始可行流 x (如零流);对节点t 标号,即令 )(fmax t =任意正值(如 1)。 STEP1 若 0)(fmax >t ,继续下一步;否则停止,已经得到最大流,结束。 STEP2 取消所有节点 Vj ∈ 的标号,即令 0)(fmax =j , 0)(pred =j ;令 LIST={ s },对节点 s 标号,即令 =)(fmax s 充分大的正值。 STEP3 如果 LIST Φ≠ 且 0)(fmax =t ,继续下一步;否则:(3a)如果t 已经有 标号(即 0)(fmax >t ),则找到了一条增广路,沿该增广路对流 x 进行增广(增广的 流量为 )(fmax t ,增广路可以根据 pred 回溯方便地得到),转 STEP1。 (3b)如果t 没有标号(即 LIST= Φ 且 0)(fmax =t ),则停止,已得到最大流。 STEP4 从 LIST 中移走一个节点i ;寻找从节点i 出发的所有可能的增广弧:(4a) 对非饱和前向弧 ),( ji ,若节点 j 没有标号(即 0)(pred =j ),对 j 进行标号,即令 },)f(min{max)(fmax ijij xuij −= , ij =)(pred , 并将 j 加入 LIST 中。 (4b)对非空后向弧 ),( ij ,若节点 j 没有标号(即 0)(pred =j ),对 j 进行标号, 即令 },)f(min{max)(fmax ijxij = , ij −=)(pred , -93- 并将 j 加入 LIST 中。 例 17 用 Ford-Fulkerson 算法计算如图 6 网络中的最大流,每条弧上的两个数字分 别表示容量和当前流量。 图 6 最大流问题 解 编写程序如下: clc,clear u(1,2)=1;u(1,3)=1;u(1,4)=2;u(2,3)=1;u(2,5)=2; u(3,5)=1;u(4,3)=3;u(4,5)=3; f(1,2)=1;f(1,3)=0;f(1,4)=1;f(2,3)=0;f(2,5)=1; f(3,5)=1;f(4,3)=1;f(4,5)=0; n=length(u);list=[];maxf(n)=1; while maxf(n)>0 maxf=zeros(1,n);pred=zeros(1,n); list=1;record=list;maxf(1)=inf; %list是未检查邻接点的标号点,record是已标号点 while (~isempty(list))&(maxf(n)==0) flag=list(1);list(1)=[]; label1= find(u(flag,:)-f(flag,:)); label1=setdiff(label1,record); list=union(list,label1); pred(label1)=flag; maxf(label1)=min(maxf(flag),u(flag,label1)... -f(flag,label1)); record=union(record,label1); label2=find(f(:,flag)); label2=label2'; label2=setdiff(label2,record); list=union(list,label2); pred(label2)=-flag; maxf(label2)=min(maxf(flag),f(label2,flag)); record=union(record,label2); end if maxf(n)>0 v2=n; v1=pred(v2); while v2~=1 if v1>0 f(v1,v2)=f(v1,v2)+maxf(n); else v1=abs(v1); f(v2,v1)=f(v2,v1)-maxf(n); end v2=v1; v1=pred(v2); end end end f -94- 例18 现需要将城市 s 的石油通过管道运送到城市t ,中间有4个中转站 321 ,, vvv 和 4v ,城市与中转站的连接以及管道的容量如图7所示,求从城市 s 到城市t 的最大流。 图7 最大流问题 解 使用最大流的数学规划表达式,编写LINGO程序如下: model: sets: nodes/s,1,2,3,4,t/; arcs(nodes,nodes)/s 1,s 3,1 2,1 3,2 3,2 t,3 4,4 2,4 t/:c,f; endsets data: c=8 7 9 5 2 5 9 6 10; enddata n=@size(nodes); !顶点的个数; max=flow; @for(nodes(i)|i #ne#1 #and# i #ne# n: @sum(arcs(i,j):f(i,j))=@sum(arcs(j,i):f(j,i))); @sum(arcs(i,j)|i #eq# 1:f(i,j))=flow; @sum(arcs(i,j)|j #eq# n:f(i,j))=flow; @for(arcs:@bnd(0,f,c)); end 在上面的程序中,采用了稀疏集的编写方法。下面介绍的程序编写方法是利用赋权邻 接矩阵,这样可以不使用稀疏集的编写方法,更便于推广到复杂网络。 model: sets: nodes/s,1,2,3,4,t/; arcs(nodes,nodes):c,f; endsets data: c=0; @text('fdata.txt')=f; enddata calc: c(1,2)=8;c(1,4)=7; c(2,3)=9;c(2,4)=5; c(3,4)=2;c(3,6)=5; c(4,5)=9;c(5,3)=6;c(5,6)=10; endcalc n=@size(nodes); !顶点的个数; max=flow; @for(nodes(i)|i #ne#1 #and# i #ne# n: @sum(nodes(j):f(i,j))=@sum(nodes(j):f(j,i))); @sum(nodes(i):f(1,i))=flow; @sum(nodes(i):f(i,n))=flow; @for(arcs:@bnd(0,f,c)); end -95- §8 最小费用流及其求法 8.1 最小费用流 上面我们介绍了一个网络上最短路以及最大流的算法,但是还没有考虑到网络上流 的费用问题,在许多实际问题中,费用的因素很重要。例如,在运输问题中,人们总是 希望在完成运输任务的同时,寻求一个使总的运输费用最小的运输方案。这就是下面要 介绍的最小费用流问题。 在运输网络 ),,,,( UAVtsN = 中,设 ijc 是定义在 A 上的非负函数,它表示通过弧 ),( ji 单位流的费用。所谓最小费用流问题就是从发点到收点怎样以最小费用输送一已 知量为 )( fv 的总流量。 最小费用流问题可以用如下的线性规划问题描述: ∑ ∈Aji ijij fc ),( min s.t. ∑∑ ∈∈ =− Aijj iji Ajij ij dff ),(:),(: , Ajiuf ijij ∈∀≤≤ ),(,0 . 其中 ⎪⎩ ⎪⎨ ⎧ ≠ =− = = tsi tifv sifv di ,,0 ),( ),( 显然,如果 =)( fv 最大流 )( maxfv ,则本问题就是最小费用最大流问题。如果 )()( maxfvfv > ,则本问题无解。 例 19(最小费用最大流问题)(续例 18)由于输油管道的长短不一或地质等原因, 使每条管道上运输费用也不相同,因此,除考虑输油管道的最大流外,还需要考虑输油 管道输送最大流的最小费用。图 8 所示是带有运费的网络,其中第 1 个数字是网络的容 量,第 2 个数字是网络的单位运费。 图 8 最小费用最大流问题 解 按照最小费用流的数学规划写出相应的 LINGO 程序如下: model: sets: nodes/s,1,2,3,4,t/:d; arcs(nodes,nodes)/s 1,s 3,1 2,1 3,2 3,2 t,3 4,4 2,4 t/:c,u,f; endsets data: d=14 0 0 0 0 -14; !最大流为14; c=2 8 2 5 1 6 3 4 7; u=8 7 9 5 2 5 9 6 10; enddata min=@sum(arcs:c*f); -96- @for(nodes(i):@sum(arcs(i,j):f(i,j))-@sum(arcs(j,i):f(j,i))=d(i)); @for(arcs:@bnd(0,f,u)); end 求得最大流的最小费用是 205,而原最大流的费用为 210 单位,原方案并不是最优 的。 类似地,可以利用赋权邻接矩阵编程求得最小费用最大流。LINGO 程序如下: model: sets: nodes/s,1,2,3,4,t/:d; arcs(nodes,nodes):c,u,f; endsets data: d=14 0 0 0 0 -14; c=0; u=0; enddata calc: c(1,2)=2;c(1,4)=8; c(2,3)=2;c(2,4)=5; c(3,4)=1;c(3,6)=6; c(4,5)=3;c(5,3)=4;c(5,6)=7; u(1,2)=8;u(1,4)=7; u(2,3)=9;u(2,4)=5; u(3,4)=2;u(3,6)=5; u(4,5)=9;u(5,3)=6;u(5,6)=10; endcalc min=@sum(arcs:c*f); @for(nodes(i):@sum(nodes(j):f(i,j))-@sum(nodes(j):f(j,i))=d(i)); @for(arcs:@bnd(0,f,u)); end 8.2 求最小费用流的一种方法—迭代法 这里所介绍的求最小费用流的方法叫做迭代法。这个方法是由 Busacker 和 Gowan 在 1961 年提出的。其主要步骤如下: (i)求出从发点到收点的最小费用通路 ),( tsμ 。 (ii)对该通路 ),( tsμ 分配最大可能的流量: }{min ),(),( ijtsji uf μ∈ = 并让通路上的所有边的容量相应减少 f 。这时,对于通路上的饱和边,其单位流费用 相应改为 ∞ 。 (iii)作该通路 ),( tsμ 上所有边 ),( ji 的反向边 ),( ij 。令 fu ji = , ijji cc −= (iv)在这样构成的新网络中,重复上述步骤(i),(ii),(iii),直到从发点到收点的全 部流量等于 )( fv 为止(或者再也找不到从 s 到t 的最小费用道路)。 下面我们编写了最小费用最大流函数 mincostmaxflow,其中调用了利用 Floyd 算法 求最短路的函数 floydpath。 求解例 19 具体程序如下(下面的全部程序放在一个文件中): function mainexample19 clear;clc; -97- global M num c=zeros(6);u=zeros(6); c(1,2)=2;c(1,4)=8;c(2,3)=2;c(2,4)=5; c(3,4)=1;c(3,6)=6;c(4,5)=3;c(5,3)=4;c(5,6)=7; u(1,2)=8;u(1,4)=7;u(2,3)=9;u(2,4)=5; u(3,4)=2;u(3,6)=5;u(4,5)=9;u(5,3)=6;u(5,6)=10; num=size(u,1);M=sum(sum(u))*num^2; [f,val]=mincostmaxflow(u,c) %求最短路径函数 function path=floydpath(w); global M num w=w+((w==0)-eye(num))*M; p=zeros(num); for k=1:num for i=1:num for j=1:num if w(i,j)>w(i,k)+w(k,j) w(i,j)=w(i,k)+w(k,j); p(i,j)=k; end end end end if w(1,num) ==M path=[]; else path=zeros(num); s=1;t=num;m=p(s,t); while ~isempty(m) if m(1) s=[s,m(1)];t=[t,t(1)];t(1)=m(1); m(1)=[];m=[p(s(1),t(1)),m,p(s(end),t(end))]; else path(s(1),t(1))=1;s(1)=[];m(1)=[];t(1)=[]; end end end %最小费用最大流函数 function [flow,val]=mincostmaxflow(rongliang,cost,flowvalue); %第一个参数:容量矩阵;第二个参数:费用矩阵; %前两个参数必须在不通路处置零 %第三个参数:指定容量值(可以不写,表示求最小费用最大流) %返回值 flow 为可行流矩阵,val 为最小费用值 global M flow=zeros(size(rongliang));allflow=sum(flow(1,:)); if nargin<3 flowvalue=M; end -98- while allflow0).*cost)'; path=floydpath(w);%调用 floydpath 函数 if isempty(path) val=sum(sum(flow.*cost)); return; end theta=min(min(path.*(rongliang-flow)+(path.*(rongliang-flow)==0).*M)); theta=min([min(path'.*flow+(path'.*flow==0).*M),theta]); flow=flow+(rongliang>0).*(path-path').*theta; allflow=sum(flow(1,:)); end val=sum(sum(flow.*cost)); §9 计划评审方法和关键路线法 计划评审方法(program evaluation and review technique, PERT)和关键路线法 (critical path method, CPM)是网络分析的重要组成部分,它广泛地用于系统分析和项 目管理。计划评审与关键路线方法是在 20 世纪 50 年代提出并发展起来的,1956 年, 美国杜邦公司为了协调企业不同业务部门的系统规划,提出了关键路线法。1958 年, 美国海军武装部在研制“北极星”导弹计划时,由于导弹的研制系统过于庞大、复杂, 为找到一种有效的管理方法,设计了计划评审方法。由于 PERT 与 CPM 既有着相同的 目标应用,又有很多相同的术语,这两种方法已合并为一种方法,在国外称为 PERT/CPM,在国内称为统筹方法(scheduling method)。 9.1 计划网络图 例 20 某项目工程由 11 项作业组成(分别用代号 KJBA ,,,, L 表示),其计划完 成时间及作业间相互关系如表 9 所示,求完成该项目的最短时间。 表 9 作业流程数据 作业 计划完成时间(天) 紧前作业 作业 计划完成时间(天) 紧前作业 A 5 - B 10 - C 11 - D 4 B E 4 A F 15 DC, G 21 EB, H 35 EB, I 25 EB, J 15 IGF ,, K 20 GF, 例 20 就是计划评审方法或关键路线法需要解决的问题。 9.1.1 计划网络图的概念 定义 称任何消耗时间或资源的行动称为作业。称作业的开始或结束为事件,事 件本身不消耗资源。 在计划网络图中通常用圆圈表示事件,用箭线表示工作,如图 9 所示,1,2,3 表 示事件, BA, 表示作业。由这种方法画出的网络图称为计划网络图。 图 9 计划网络图的基本画法 -99- 虚工作用虚箭线“……→”表示。它表示工时为零,不消耗任何资源的虚构工作。 其作用只是为了正确表示工作的前行后继关系。 定义 在计划网络图中,称从初始事件到最终事件的由各项工作连贯组成的一条 路为路线。具有累计作业时间最长的路线称为关键路线。 由此看来,例 20 就是求相应的计划网络图中的关键路线。 9.1.2 建立计划网络图应注意的问题 (1)任何作业在网络中用唯一的箭线表示,任何作业其终点事件的编号必须大于 其起点事件。 (2)两个事件之间只能画一条箭线,表示一项作业。对于具有相同开始和结束事 件的两项以上的作业,要引进虚事件和虚作业。 (3)任何计划网络图应有唯一的最初事件和唯一的最终事件。 (4)计划网络图不允许出现回路。 (5)计划网络图的画法一般是从左到右,从上到下,尽量作到清晰美观,避免箭 头交叉。 9.2 时间参数 9.2.1 事件时间参数 (1)事件的最早时间 事件 j 的最早时间用 )( jtE 表示,它表明以它为始点的各工作最早可能开始的时 间,也表示以它为终点的各工作的最早可能完成时间,它等于从始点事件到该事件的最 长路线上所有工作的工时总和。事件最早时间可用下列递推公式,按照事件编号从小到 大的顺序逐个计算。 设事件编号为 n,,2,1 L ,则 ⎪⎩ ⎪⎨ ⎧ += = )},()({max)( 0)1( jititjt t EiE E (6) 其中 )(itE 是与事件 j 相邻的各紧前事件的最早时间, ),( jit 是作业 ),( ji 所需的工时。 终点事件的最早时间显然就是整个工程的总最早完工期,即 =)(ntE 总最早完工期 (7) (2)事件的最迟时间 事件i 的最迟时间用 )(itL 表示,它表明在不影响任务总工期条件下,以它为始点 的工作的最迟必须开始时间,或以它为终点的各工作的最迟必须完成时间。由于一般情 况下,我们都把任务的最早完工时间作为任务的总工期,所以事件最迟时间的计算公式 为: ⎪⎩ ⎪⎨ ⎧ −= = )},()({min)( )()( jitjtit ntnt LjL EL )总工期(或 (8) 其中 )( jtL 是与事件i 相邻的各紧后事件的最迟时间。 公式(8)也是递推公式,但与(6)相反,是从终点事件开始,按编号由大至小的 顺序逐个由后向前计算。 9.2.2 工作的时间参数 (1)工作的最早可能开工时间与工作的最早可能完工时间 一个工作 ),( ji 的最早可能开工时间用 ),( jitES 表示。任何一件工作都必须在其所 -100- 有紧前工作全部完工后才能开始。工作 ),( ji 的最早可能完工时间用 ),( jitEF 表示。它 表示工作按最早开工时间开始所能达到的完工时间。它们的计算公式为: ⎪⎪ ⎩ ⎪⎪ ⎨ ⎧ += += = ),(),(),( )},(),({max),( 0),1( jitjitjit iktiktjit jt ESEF ESkES ES (9) 这组公式也是递推公式。即所有从总开工事件出发的工作 ),1( j ,其最早可能开工 时间为零;任一工作 ),( ji 的最早开工时间要由它的所有紧前工作 ),( ik 的最早开工时间 决定;工作 ),( ji 的最早完工时间显然等于其最早开工时间与工时之和。 (2)工作的最迟必须开工时间与工作的最迟必须完工时间 一个工作 ),( ji 的最迟开工时间用 ),( jitLS 表示。它表示工作 ),( ji 在不影响整个任 务如期完成的前提下,必须开始的最晚时间。 工作 ),( ji 的最迟必须完工时间用 ),( jitLF 表示。它表示工作 ),( ji 按最迟时间开 工,所能达到的完工时间。它们的计算公式为: ⎪⎪ ⎩ ⎪⎪ ⎨ ⎧ += −= = ),(),(),( )},(),({min),( ),(),( jitjitjit jitkjtjit nitnit LSLF LSkLS EFLF )总完工期(或 (10) 这组公式是按工作的最迟必须开工时间由终点向始点逐个递推的公式。凡是进入 总完工事件 n 的工作 ),( ni ,其最迟完工时间必须等于预定总工期或等于这个工作的最 早可能完工时间。任一工作 ),( ji 的最迟必须开工时间由它的所有紧后工作 ),( kj 的最 迟开工时间确定。而工作 ),( ji 的最迟完工时间显然等于本工作的最迟开工时间与工时 的和。 由于任一个事件i (除去始点事件和终点事件),既表示某些工作的开始又表示某 些工作的结束。所以从事件与工作的关系考虑,用公式(9),公式(10)求得的有关工 作的时间参数也可以通过事件的时间参数公式(6),公式(8)来计算。如工作 ),( ji 的 最早可能开工时间 ),( jitES 就等于事件i 的最早时间 )(itE 。工作 ),( ji 的最迟必须完工 时间等于事件 j 的最迟时间。 9.2.3 时差 工作的时差又叫工作的机动时间或富裕时间,常用的时差有两种。 (1)工作的总时差 在不影响任务总工期的条件下,某工作 ),( ji 可以延迟其开工时间的最大幅度,叫 做该工作的总时差,用 ),( jiR 表示。其计算公式为: ),(),(),( jitjitjiR EFLF −= (11) 即工作 ),( ji 的总时差等于它的最迟完工时间与最早完工时间的差。显然 ),( jiR 也等于 该工作的最迟开工时间与最早开工时间之差。 (2)工作的单时差 工作的单时差是指在不影响紧后工作的最早开工时间条件下,此工作可以延迟其 开工时间的最大幅度,用 ),( jir 表示。其计算公式为: -101- ),(),(),( jitkjtjir EFES −= (12) 即单时差等于其紧后工作的最早开工时间与本工作的最早完工时间之差。 9.3 计划网络图的计算 以例 20 的求解过程为例介绍计划网络图的计算方法。 9.3.1 建立计划网络图 首先建立计划网络图。按照上述规则,建立例 20 的计划网络图,如图 10 所示。 图 10 例 20 的计划网络图 9.3.2 写出相应的规划问题 设 ix 是事件i 的开始时间,1为最初事件, n 为最终事件。希望总的工期最短,即 极小化 1xxn − 。设 ijt 是作业 ),( ji 的计划时间,因此,对于事件i 与事件 j 有不等式 ijij txx +≥ 由此得到相应的数学规划问题 1min xxn − s.t. ijij txx +≥ , Aji ∈),( , Vji ∈, 0≥ix , Vi ∈ 其中V 是所有的事件集合, A 是所有的作业集合。 9.3.3 问题求解 用 LINGO 软件求解例 20。 解 编写 LINGO 程序如下: model: sets: events/1..8/:x; operate(events,events)/1 2,1 3,1 4,2 5,3 4,3 5,4 6,5 6, 5 7,5 8,6 7,6 8,7 8/:t; endsets data: t=5 10 11 4 4 0 15 21 25 35 0 20 15; enddata min=x(8)-x(1); @for(operate(i,j):x(j)>x(i)+t(i,j)); end 计算结果给出了各个项目的开工时间,如 01 =x ,则作业 CBA ,, 的开工时间均是 第 0 天; 52 =x ,作业 E 的开工时间是第 5 天; 103 =x ,则作业 D 的开工时间是第 10 天;等等。每个作业只要按规定的时间开工,整个项目的最短工期为 51 天。 尽管上述 LINGO 程序给出相应的开工时间和整个项目的最短工期,但统筹方法中 -102- 许多有用的信息并没有得到,如项目的关键路径、每个作业的最早开工时间、最迟开工 时间等。 例 21(续例 20)求例 20 中每个作业的最早开工时间、最迟开工时间和作业的关 键路径。 解 为了得到每个作业的最早开工时间、作业的关键路线等,将目标函数改为 ∑ ∈Vi ix ,即作业的开始时间尽量早,这样就可以得到作业的最早开工时间。再引进作业 对应弧上的松弛变量 ijs ,且 ijijij txxs −−= , Aji ∈),( ,这样就可以得到作业的最迟 开工时间,记 iy 表示事件i 的最迟开工时间。当最早开工时间与最迟开工时间相同时, 就得到项目的关键路径。 编写 LINGO 程序如下: model: sets: events/1..8/:x,z; operate(events,events)/1 2,1 3,1 4,2 5,3 4,3 5,4 6,5 6, 5 7,5 8,6 7,6 8,7 8/:s,t,m,c,y; endsets data: t=5 10 11 4 4 0 15 21 25 35 0 20 15; m=5 8 8 3 4 0 15 16 22 30 0 16 12; c=0 700 400 450 0 0 0 600 300 500 0 500 400; d=49; @text(txt2.txt)=x,z; enddata min=mincost+sumx; mincost=@sum(operate:c*y); sumx=@sum(events:x); @for(operate(i,j):s(i,j)=x(j)-x(i)+y(i,j)-t(i,j)); n=@size(events); x(1)=0; x(n)ijs 时,说明还有剩余时间,对 应作业的工期可以推迟 ijs 。例如, 178 =s ,作业(7,8)( J )的开工时间可以推迟 1 天,即开工时间为 36。再如 246 =s ,作业(4,6)( F )可以推迟 2 天开始, 314 =s , 作业(1,4)(C )可以推迟 3 天开始,但由于作业(4,6)( F )已能够推迟 2 天, 所以,作业(1,4)(C )最多可推迟 5 天。 由此,可以得到所有作业的最早开工时间和最迟开工时间,如下表所示,方括号 中第 1 个数字是最早开工时间,第 2 个数字是最迟开工时间。 表 10 作业数据 作业 ),( ji 开工时间 计划完成时间(天) 作业 ),( ji 开工时间 计划完成时间(天) )2,1(A [0,1] 5 )6,5(G [10,10] 21 -103- )3,1(B [0,0] 10 )4,1(C [0,5] 11 )4,3(D [10,12] 4 )5,2(E [5,6] 4 )6,4(F [14,16] 15 )8,5(H [10,16] 35 )7,5(I [10,11] 25 )8,7(J [35,36] 15 )8,6(K [31,31] 20 从上表可以看出,当最早开工时间与最迟开工时间相同时,对应的作业在关键路 线上,因此可以画出计划网络图中的关键路线,如图 11 粗线所示。关键路线为 1→3→ 5→6→8。 图 11 带有关键路线的计划网络图 9.3.4 将关键路线看成最长路 如果将关键路线看成最长路,则可以按照求最短路的方法(将求极小改为求极大) 求出关键路线。 设 ijx 为 0-1 变量,当作业 ),( ji 位于关键路线上取 1,否则取 0。数学规划问题写 成: ∑ ∈Aji ijij xt ),( max s.t. ⎪⎩ ⎪⎨ ⎧ ≠ =− = =− ∑∑ ∈ = ∈ = ni ni i xx n Aij j ji n Aji j ij ,1,0 ,1 1,1 ),( 1 ),( 1 0=ijx 或 1, Aji ∈),( 例 22 用最长路的方法,求解例 20。 解 按上述数学规划问题写出相应的 LINGO 程序。 model: sets: events/1..8/:d; operate(events,events)/1 2,1 3,1 4,2 5,3 4,3 5,4 6,5 6, 5 7,5 8,6 7,6 8,7 8/:t,x; endsets data: t=5 10 11 4 4 0 15 21 25 35 0 20 15; d=1 0 0 0 0 0 0 -1; enddata max=@sum(operate:t*x); @for(events(i):@sum(operate(i,j):x(i,j))-@sum(operate(j,i):x(j,i))=d(i)); end -104- 求得工期需要 51 天,关键路线为 1→3→5→6→8。 9.4 关键路线与计划网络的优化 例 23(关键路线与计划网络的优化) 假设例 20 中所列的工程要求在 49 天内完成。 为提前完成工程,有些作业需要加快进度,缩短工期,而加快进度需要额外增加费用。 下表列出例 20 中可缩短工期的所有作业和缩短一天工期额外增加的费用。现在的问题 是,如何安排作业才能使额外增加的总费用最少。 表 11 工程作业数据 作业 计划完成 最短完成 缩短 1 天增加 ),( ji 时间(天) 时间(天) 的费用(元) 作业 计划完成 最短完成 缩短 1 天增加 ),( ji 时间(天) 时间(天) 的费用(元) )3,1(B 10 8 700 )4,1(C 11 8 400 )5,2(E 4 3 450 )6,5(G 21 16 600 )8,5(H 35 30 500 )7,5(I 25 22 300 )8,7(J 15 12 400 )8,6(K 20 16 500 例 23 所涉及的问题就是计划网络的优化问题,这时需要压缩关键路径来减少最短 工期。 9.4.1 计划网络优化的数学表达式 设 ix 是事件i 的开始时间, ijt 是作业 ),( ji 的计划时间, ijm 是完成作业 ),( ji 的最 短时间, ijy 是作业 ),( ji 可能减少的时间, ijc 是作业 ),( ji 缩短一天增加的费用,因此 有 ijijij ytxx −≥− 且 ijijij mty −≤≤0 设 d 是要求完成的天数,1 为最初事件,n 为最终事件,所以有 dxxn ≤− 1 。而问题的 总目标是使额外增加的费用最小,即目标函数为 ∑ ∈Aji ijij yc ),( min 。由此得到相应的数学 规划问题: ∑ ∈Aji ijij yc ),( min s.t. ijijij tyxx ≥+− , Aji ∈),( , Vji ∈, dxxn ≤− 1 ijijij mty −≤≤0 , Aji ∈),( , Vji ∈, 9.4.2 计划网络优化的求解 用 LINGO 软件求解例 23,程序如下: model: sets: events/1..8/:x; operate(events,events)/1 2,1 3,1 4,2 5,3 4,3 5,4 6,5 6, 5 7,5 8,6 7,6 8,7 8/:t,m,c,y; endsets data: t=5 10 11 4 4 0 15 21 25 35 0 20 15; -105- m=5 8 8 3 4 0 15 16 22 30 0 16 12; c=0 700 400 450 0 0 0 600 300 500 0 500 400; d=49; enddata min=@sum(operate:c*y); @for(operate(i,j):x(j)-x(i)+y(i,j)>t(i,j)); n=@size(events); x(n)-x(1)= 500*f(i)); @for(need(j):[con3] @sum(supply(i):x(i,j))=y(j)+z(j)); @for(need(j)|j#NE#15:[con4] z(j)+y(j+1)=b(j)); y(1)=0; z(15)=0; @for(supply: @bin(f)); @for(need: @gin(y)); end 10.3 管道为树形图时的模型 当管道为树形图时,建立与上面类似的非线性规划模型 )(05.0min 2 21 1)( 7 1 21 1 jkjk jEjk ij ij ij yyxc ++ ∑∑∑∑ =∈== (27) s.t. ∑ = ≤≤ 21 1 500 j iiiji fsxf , 7,,2,1 L=i (28) ∑∑ =∈ = 7 1)(i jk