多目标规划MAABwg.ppt

MATLAB求解多目标规划 江西师范大学数信学院吴根秀 一 0 1规划的MATLAB求解 数学模型 MINf xS T Ax bAeqx beqx 0 1命令格式 x bintprog f x bintprog f A b x bintprog f A b Aeq beq x bintprog f A b Aeq beq x0 x bintprog f A b Aeq beq x0 options x fval bintprog x fval exitflag bintprog x fval exitflag output bintprog 数学模型 MINlambdaS T F x weight lambda goal 达到目标 Ax b 线性不等式约束 Aeqx beq 线性等式约束 C x 0 非线性不等式约束 Ceq x 0 非线性等式约束 lb x ubF f1 x f2 x 为多目标的目标函数 F与 C x Ceq x 都是通过function来定义 命令格式 x fgoalattain fun x0 goal weight x fgoalattain fun x0 goal weight A b x fgoalattain fun x0 goal weight A b Aeq beq x fgoalattain fun x0 goal weight A b Aeq beq lb ub 二 多目标规划的MATLAB求解 命令格式 x fgoalattain fun x0 goal weight x fgoalattain fun x0 goal weight A b x fgoalattain fun x0 goal weight A b Aeq beq x fgoalattain fun x0 goal weight A b Aeq beq lb ub x fgoalattain fun x0 goal weight A b Aeq beq lb ub nonlcon x fgoalattain fun x0 goal weight A b Aeq beq lb ub nonlcon options x fval fgoalattain x fval attainfactor fgoalattain x fval attainfactor exitflag fgoalattain x fval attainfactor exitflag output fgoalattain x fval attainfactor exitflag output lambda fgoalattain 二 多目标规划的MATLAB求解 x fgoalattain myfun x0 goal weight A b Aeq beq lb ub mycon wheremyconisaMATLABfunctionsuchasfunction c ceq mycon x c computenonlinearinequalitiesatx ceq computenonlinearequalitiesat 二 多目标规划的MATLAB求解 x fgoalattain myfun x0 goal weight wheremyfunisaMATLABfunctionsuchasfunctionF myfun x F Computefunctionvaluesatx 有关优化参数设置 options optimset GradObj on 目标函数的梯度方向参数设置为 on 时 用下列函数定义 function F G myfun x F Computethefunctionvaluesatxifnargout 1 TwooutputargumentsG GradientsevaluatedatxEndThegradientconsistsofthepartialderivativedF dxofeachFatthepointx 二 多目标规划的MATLAB求解 二 多目标规划的MATLAB求解 有关优化参数设置 options optimset GradConstr on 约束条件的梯度方向参数设置为 on 时 用下列函数定义 function c ceq GC GCeq mycon x c Nonlinearinequalitiesatxceq Nonlinearequalitiesatxifnargout 2 Nonlconcalledwith4outputsGC GradientsoftheinequalitiesGCeq GradientsoftheequalitiesEnd注意 一般weight abs goal 模型 x A BKC x Bu 设计K满足目标 Y Cx1 循环系统的特征值 由命令eig A B K C 确定 的目标为goal 5 3 1 2 K中元素均在 4 4 中 设特征值的weight abs goal 定义目标函数F如下 functionF eigfun K A B C F sort eig A B K C Evaluateobjectives 由小到大排列优化程序为 A 0 500 0 210 01 2 B 10 22 01 C 100 001 K0 1 1 1 1 Initializecontrollermatrixgoal 5 3 1 Setgoalvaluesfortheeigenvaluesweight abs goal Setweightforsamepercentagelb 4 ones size K0 Setlowerboundsonthecontrollerub 4 ones size K0 Setupperboundsonthecontrolleroptions optimset Display iter Setdisplayparameter K fval attainfactor fgoalattain K eigfun K A B C goal weight lb ub options 二 举例 有关循环控制系统优化问题 运行结果如下Activeconstraints 124910K 4 0000 0 2564 4 0000 4 0000fval 6 9313 4 1588 1 4099attainfactor 0 3863 二 举例 有关循环控制系统优化问题 如果至少保证38 63 的目标精确匹配 设置 GoalsExactAchieve 参数值为3options optimset GoalsExactAchieve 3 K fval attainfactor fgoalattain K eigfun K A B C K0 goal weight lb ub options Afteraboutseveniterations asolutionisK 1 59541 2040 0 4201 2 9046fval 5 0000 3 0000 1 0000attainfactor 1 0859e 20表明目标已完全匹配 二 举例 有关循环控制系统优化问题 谢谢 初等模型举例 常见类型 定性模型经验公式 拟合 插值 量纲分析比例模型 2 1崖高的估算 假如你站在崖顶且身上带着一只具有跑表功能的计算器 你也许会出于好奇心想用扔下一块石头听回声的方法来估计山崖的高度 假定你能准确地测定时间 你又怎样来推算山崖的高度呢 请你分析一下这一问题 方法一 我学过微积分 我可以做得更好 呵呵 令k K m 解得 代入初始条件v 0 0 得c g k 故有 再积分一次 得 若设k 0 05并仍设t 4秒 则可求得h 73 6米 听到回声再按跑表 计算得到的时间中包含了反应时间 进一步深入考虑 不妨设平均反应时间为0 1秒 假如仍设t 4秒 扣除反应时间后应为3 9秒 代入式 求得h 69 9米 多测几次 取平均值 再一步深入考虑 2 2录像带还能录多长时间 录像机上有一个四位计数器 一盘180分钟的录像带在开始计数时为0000 到结束时计数为1849 实际走时为185分20秒 我们从0084观察到0147共用时间3分21秒 若录像机目前的计数为1428 问是否还能录下一个60分钟的节目 又因和得 积分得到 即 从而有 此式中的三个参数W v和r均不易精确测得 虽然我们可以从上式解出t与n的函数关系 但效果不佳 故令则可将上式简化为 故 t an2 bn 上式以a b为参数显然是一个十分明智的做法 它为公式的最终确立即参数求解提供了方便 将已知条件代入 得方程组 从后两式中消去t1 解得a 0 0000291 b 0 04646 故t 0 0000291n2 0 04646n 令n 1428 得到t 125 69 分 由于一盒录像带实际可录像时间为185 33分 故尚可录像时间为59 64分 已不能再录下一个60分钟的节目了 2 3最短路径与最速方案问题 例5 最短路径问题 设有一个半径为r的圆形湖 圆心为O A B位于湖的两侧 AB连线过O 见图 现拟从A点步行到B点 在不得进入湖中的限制下 问怎样的路径最近 以上只是一种猜测 现在来证明这一猜测是正确的 为此 先介绍一下凸集与凸集的性质 下面证明猜想 猜测证明如下 还可用微积分方法求弧长 根据计算证明满足限止条件的其他连续曲线必具有更大的长度 此外 本猜测也可用平面几何知识加以证明等 到此为止 我们的研讨还只局限于平面之中 其实上述猜测可十分自然地推广到一般空间中去 1973年 J W Craggs证明了以上结果 例6一辆汽车停于A处并垂直于AB方向 此汽车可转的最小圆半径为R 求不倒车而由A到B的最短路径 例7驾驶一辆停于A处与AB成 1角度的汽车到B处去 已知B处要求的停车方向必须与AB成 2角 试找出最短路径 除可转的最小圆半径为R外 不受其他限止 最速方案问题 例8将一辆急待修理的汽车由静止开始沿一直线方向推至相隔S米的修车处 设阻力不计 推车人能使车得到的推力f满足 B f A f 0为推力 f 0为拉力 问怎样推车可使车最快停于修车处 此问题为一泛函极值问题 求解十分困难 为得出一个最速方案 我们作如下猜测 逻辑模型 例1拟将一批尺寸为1 2 4的的商品装入尺寸为6 6 6的正方体包装箱中 问是否存在一种装法 使装入的该商品正好充满包装箱 解将正方体剖分成27个2 2 2的小正方体 并按下图所示黑白相间地染色 再将每一2 2 2的小正方体剖分成1 1 1的小正方体 易见 27个2 2 2的正方体中 有14个是黑的 13个是白的 或13黑14白 故经两次剖分 共计有112个1 1 1的黑色小正方体和104个1 1 1的白色小正方体 虽然包装箱的体积恰好是商品体积的27倍 但容易看到 不论将商品放置在何处 它都将占据4个黑色和4个白色的1 1 1小正方体的位置 故商品不可能充满包装箱 德国著名的艺术家AlbrechtD rer 1471 1521 于1514年曾铸造了一枚名为 MelencotiaI 的铜币 令人奇怪的是在这枚铜币的画面上充满了数学符号 数字及几何图形 这里 我们仅研究铜币右上角的数字问题 所谓的魔方是指由1 n2这n2个正整数按一定规则排列成的一个n行n列的正方形 n称为此魔方的阶 D rer魔方 4阶 每一行之和为34 每一列之和为34 对角线 或反对角线 之和是34 每个小方块中的数字之和是34 四个角上的数字加起来也是34 什么是D rer魔方 多么奇妙的魔方 铜币铸造时间 1514年 构造魔方是一个古老的数学游戏 起初它还和神灵联系在一起 带有深厚的迷信色彩 传说三千二百多年前 公元前2200年 因治水出名皇帝大禹就构造了三阶魔方 被人们称 洛书 至今还有人把它当作符咒用于某些迷信活动 大约在十五世纪时 魔方传到了西方 著名的科尼利厄斯 阿格里帕 1486 1535 先后构造出了3 9阶的魔方 如何构造魔方 奇数 不妨n 5 阶的情况 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 魔方数量随阶数n增长的速度实在是太惊人了 同阶魔方的个数 允许构成魔方的数取任意实数 允许取实数 n阶魔方A B 任意实数 A B是n阶魔方 具有指定性质的魔方全体构成一个线性空间 问题已发生了实质性变化 注 刻画一个线性空间只需指出它的维数并求出此线性空间的一组基底 松驰问题的讨论 1在第一行中共有4种取法 为保持上述性质的成立 第二行中的1还有两种取法 当第二行的1也取定后 第三行与第四行的1就完全定位了 故一共可作出8个不同的最简方阵 称之为基本魔方并记之为Q1 Q8 仍以4阶方阵为例 令R为行和 C为列和 D为对角线和 S为小方块和定义0 方 R C D S 0定义1 方 R C D S 4 R C D S 1的方阵构成的线性空间具有什么样的性质 类似于构造n维欧氏空间的标准基 利用0和1我们来构造一些R C D S 1的最简单的方阵 显然 D rer空间 简称D空间 中任何一个元素都可以用Q1 Q2 Q8来线性表示 但它们能否构成D空间的一组基呢 等号两边对应元素相比较 得r1 r2 r7 0 所以是线性无关 是D空间的最小生成集 令D即解方程组 解得D 研究AlbrechtD rer铸造的铜币 2004年浙江大学数学建模竞赛 B题 通讯卫星上的开关设置地面上存在着n个接收站与n个发送站 而在通讯卫星上则设置了若干种开关模式 开关模式可用矩阵P pij 来表示 若卫星可接收发送站i发射的信息并将信息传送回地面的接收站j时 矩阵中的元素pij 1 否则pij 0 通讯卫星上的接收发送任务也可以用一个矩阵T tij 来表示 其元素tij为需经通讯卫星传递的由i发点发送到j接收点的信息量的传送时间长度 由于技术上的原因 当发送站i在发送给接收站j信息时 它不能同时发送给别的接收站信息 同样 当接收站j在接收发送站i的信息时 也不能同时接收其他发送站发送的信息 你的任务是 设计一组开关模式 k 1 r 注 r应当尽可能小 使得对任意给定的任务矩阵T 卫星开关设置均能完成要求的发接收任务 设计一个算法 在发接收任务T给出后 可根据你设计的开关模式 k 1 r 求出各开关的使用时间 k 使得在完成预定传送任务的前提下使用各开关模式的总时间最短 同样由于技术上的原因 开关模式的总数r有一个上限 当需要传送的任务数数量较大时 仍无法分派任务 请你想一些办法来解决这一困难 当然 这时你可能要作出一些牺牲 即传送时间可能会增加一些 问题及模型 问题的标准形式为 在地面上存在着n个收站与n个发战 而在通讯卫星上则设置了若干种开关模式 开关模式可用矩阵P pij 来表示 若卫星可接收发射站i发射的信息并将信息传送回地面的接收站j时 矩阵元素pij 1 否则pij 0 通讯卫星的接发任务也可用一矩阵T tij 来表示 其元素tij为需经通讯卫星传递的由i发点发送到j接受点的信息量的传送时间长度 问题要求求r并设计一组开关模式Pk k 1 r及模式Pk的使用时间 k 使得在完成预定传送任务的前提下各开关模式使用的总时间最短 即要求求解下面的问题 例1设 这是一个有3个发送站与3个接收站的实例 tij在矩阵中已给出 例如由发站1传送到收站1的通讯量为3单位时间等 分析容易看出 三个发站需传送的时间分别为6 5 5 而三个收站需接收的时间分别为6 3 7 为完成全部传送任务 通讯卫星总传送时间至少应为7单位时间 即的下界为7 由于技术上的原因 当发站i在发送给收站j信息时 它不能同时发送给别的收站信息 同样 当收站j在接收发站i的信息时 也不能同时接收其他发站发送的信息 这一要求说明 任一开关模式Pk应具有以下性质 1 Pk的每一行中有且只有一个1 每一列中也有且只有一个1 2 所有的1均位于不同的行列中 满足 1 2 的矩阵被称为置换矩阵 n阶置换矩阵Pk共有n 个 当n较大时 我们不可能在通讯卫星上设置这么多种不同的开关模式 因而 为了设计出切实可行的开关模式 我们还得另想办法 设计方法1 注意到Pk每行 或列 元素之和均为1 故不管如何指派开关的使用时间 即不论如何取 k 矩阵 均具有某些特殊的性质 例如其行和 及列和 均为同一常数 这样的矩阵构成一个线性空间 参见逻辑模型第一节D rer魔方 为减少开关模式的种类 可取此空间的一组基底作为开关模式 在使用这种开关模式时 无论T的元素tij怎么取 通讯卫星对每一发 收 点的开通时间总和是恒定的 在这种开关模式下 可按如下方式指派各开关模式的使用时间 步1先将T改变为 满足 将T化为的方法一般有无穷多种 如可如下化法 令事实上 即通讯卫星传送总时间的下界 令 其中 用这种方法化例中的T 得到 的任一行 或列 中元素之和均为7 定义1称行和 列 和均相等的矩阵为双随机矩阵 Doublystochasticmatrix 定理1 Birkhoff定理 1944 任一n阶双随机矩阵均可写成至多 n 1 2 1个置换矩阵的非线性组合 的分解可如下进行 步1选取由Pij 0可推出 0的置换矩阵P 步2确定 步3取 用 代替 步4若 0 停 否则 返回步1 例2 为方便起见 我们来分解一个元素均为非负整数的3阶双随机矩阵 由Birkhoff定理 r 5 解 取 min 1 3 3 1 因min 5 5 3 3 又有 取 于是又有 易得分解结果为 尚需解决的问题是如何求P 使得Pij 0必有 读者不难发现 此问题可以通过求解一个两分图上的最大流 或最大匹配 来实现 计算量为O n4 是多项式时间可解的 具体方法为 作一两分图 若 则作边 i j 令边容量为1 这样 可作出P的充要条件是该最大流问题的最大流量为n 对例9 33 n 3 由于所有 先取 P1为 相应的两分图为 于是又可求得 相应的两分图为 又可得 如此下去 直到作不出P为至 由于的特殊性质及Birkhoff定理 上述分解必能在不超过r n 1 2 1步内终止 上述开关设计方法要求在通讯卫星上设置 n 1 2 1种不同的开关模式 即Pk 当n稍大时 n 1 2 1仍显得太大而使得使用时不便 例如 当n 41时 n 1 2 1 60 为实用方便 人们研究了限止开关模式个数的相应问题 若要求r n 即要求通讯卫星上至多设置n种开关模式 则问题化为令r n 求不超过n个置换矩阵Pk及 k 使之满足 minS t 为了使任意一对发射法与接收站之间的传送均为可能实现的 自然应要求Pk满足 5 1 5 2 右面的矩阵有n2个值为1的分量 每一Pk恰有n个1分量 故r n 容易看出 5 1 隐含着T的每一元素只能被唯一的P复盖 即T的元素在分解中是不可分割的 这当然是一个好性质 使实际操作时较为方便 但可惜的是对一般的双随机矩阵 分解很可能无解 例3若取 注意 T已是双随机矩阵 行和列和均为10 则min S t 的解为 1 3 2 4 3 5 大于10 而 但等号经常并不成立 1985年 F Rendel证明 在给定满足 5 2 的置换矩阵P1 Pn后 求解问题 5 1 是NP难的 从而不可能存在多项式时间算法 除非P NP 现要求r 2n 一种自然而方便的开关设置为引入两组各有n个开关模式的置换矩阵P1 Pn Q1 Qn 满足下面的 5 3 式 例如 当n 3时 可令 注 这种设置方法保持了其内在的对称性 不失为一种明智的做法 现在 我们来分解例9 33中的双随机矩阵 令 得方程组 求出各对角线与反对角线上的三个元素之和 并作一些简单的消去运算 将矩阵的所有元素相加 可得下面的方程组 注意到 5 3 易证空间的维数为5 故之一可任取 稍加注意即可保持非负性 例如 令 3 0 求得 故有 读者不难验证 上述方法可推广到n是奇数的一般情况 事实 由各