数学建模数学规划模型课件.ppt

数学规划模型的一般表达式： 为目标函数,为约束函数， 为约束函数， 为可控变量， 为已知参数， 为随机参数。 数学规划分为线性规划、非线性规划、动态规划、 随机规划、整数规划、分式规划、几何规划、目标规划、平衡规划、参数规划、多目标规划等十几种。当然这么多规划其中亦有交叉。又可经过组产生新的规划，每一种规划有专著问世。 min max,f x . .,0s t g x fgx第一节 线性规划模型), 2 , 1(0.min221122222121112121112211nixbxaxaxabxaxaxabxaxaxatsxcxcxcfimnmnmmnnnnnn（1）目标函数是决策变量的线性函数。 （2）约束条件都是决策变量的线性等式或不等式。 MATLAB命令命令输入格式及线性规划模型如下：其中：x0是算法迭代的初始点；nEq表示等式约束的个数。), 0,(nEqxxUBxLBbAclpX 2211.minbxAbxAtsxcfxUBxxLB2121,bbbAAA三、建模举例营养配餐问题 每种蔬菜含有的营养素成份是不同的，从医学上知道，每人每周对每种营养成分的最低需求量。某医院营养室在制定下一周菜单时，需要确定表6-1中所列六种蔬菜的供应量，以便使费用最小而又能满足营养素等其它方面的要求。规定白菜的供应一周内不多于20千克，其它蔬菜的供应在一周内不多于40千克，每周共需供应140千克蔬菜，为了使费用最小又满足营养素等其它方面的要求，问在下一周内应当供应每种蔬菜各多少千克？表2-3问题分析与模型建立设 分别表示下一周内应当供应的青豆、胡萝卜、菜花、白菜、甜菜及土豆的量，费用目标函数为：约束条件：铁的需求量至少6个单位数：磷的需求量至少25个单位数：(1 6)ix i 123456558263fxxxxxx1234560.450.451.050.400.500.506xxxxxx12345610285925227525xxxxxx维生素A的需求量至少17500个单位：维生素C的需求量至少245个单位：烟酸的需求量至少5个单位数：每周需供应140千克蔬菜，即12345641590652550751523517500 xxxxxx12345683532758245xxxxxx1234560.300.350.600.150.250.805xxxxxx123456140 xxxxxx 0 x140 0 x240 0 x340 0 x420 0 x540 0 x640123456123456123456123456123456123456min5582631400.450.451.050.400.500.50610285925227525.41590652550751523517500835327582450.3fxxxxxxxxxxxxxxxxxxxxxxxxstxxxxxxxxxxxx12345600.350.600.150.250.805xxxxxx问题是满足营养素要求的条件下，所需费用最小，是一个线性规划模型。 利用Matlab软件编程序： % 营养配餐ch21% 文件名： ch21 mc=5;5;8;2;6;3;A=(-1)*1,1,1,1,1,1; 0.45,0.45,1.05,0.40,0.50,0.50; 10,28,59,25,22,75; 415,9065,2550,75,15,235;8,3,53,27,5,8; 0.30,0.35,0.60,0.15,0.25,0.80;b=(-1)*140;6;25;17500;245;5;xLB=zeros(6,1);xUB=40;40;40;20;40;40;nEq=1;x0=0*ones(6,1);x=lp(c,A,b,xLB,xUB,x0,nEq);disp(青豆需要的份数)x(1)disp(胡罗卜需要的份数 ) x(2)disp(菜花需要的份数 )x(3)disp(白菜需要的份数 )x(4)disp(甜菜需要的份数 )x(5)disp(土豆需要的份数)x(6)执行后输出青豆需要的份数ans = 40胡罗卜需要的份数 ans = 40.0000菜花需要的份数 ans = 0白菜需要的份数 ans = 20.0000甜菜需要的份数 ans = 0土豆需要的份数ans = 40最小费用 ans = 560.0000 背景：0-1规划是数学规划的组成部分，起始20世纪30年代末，七八十年代是数学规划飞速发展时期，无论是从理论上还是算法方面都得到了进一步完善。时至今日数学规划仍然是运筹学领域中热点研究问题，从国内外的数学建模竞赛的试题中看，有近1/2的问题可用数学规划进行求解。其中利用0-1 规划及0-1型变量的数学建模问题也为数不少，如98年的投资的收益和风险 ，2004年的DVD在线租赁等问题，下面我们就来学习0-1规划, 0-1型变量在数学建模中的应用。 2.2 0-1规划, 0-1型变量 在数学建模中的应用 1、 0-1规划规划 数学规划模型的一般表达式： 整数规划中决策变量只取0或1的特殊情况是0-1规划。下面通过几个例子说明0-1规划在实际问题中的应用。例例2.1 背包问题背包问题 有几件物品，编号为 1，2，n。第 件重为 kg，价值为 元。今有一位装包者欲将这些物品装入一包，其质量不能超过 kg，问应装入哪几件价值最大？iiaipa ,(,.,)12min maxxxxxnf x其中 . .0st g x 解解 引入变量 ， 将 物品装包 ， 不将 物品装包 于是得问题的模型为 取0或1，i=1，2，n 背包问题看似简单，但应用很广，例如某些投资问题即可归入背包问题模型。此类问题可以描述为：10ixii11 1max. .niiinnp xsta xa xaix,ix,ix 投资问题：投资问题：设有总额为 元的资金，投资几项事业，第 项事业需投资 元，利润为 元，问应选择哪些项投资总利润为最大？例例2.2 某钻井队要从以下10个可供选择的井位中确定5个钻井探油，使总的钻探费用为最小。若10个井位的代号为 ，相应的钻探费用为 ，并且井位选择要满足下列限制条件： （1）或选 和 ，或选 ； （2）选择了 或 就不能选 ，反之亦然； （3）在 中最多只能选两个。 试建立其数学模型：aiiaib1210,s ss1210,c cc1s7s8s3s4s5s5678,s s s s解解 引入变量 选择 不选择 于是以上问题的数学模型为：10ixisis1 01m iniiic xix1 01187835455678511.11201iiixxxxxstxxxxxxxxx 投资的收益和风险 这是1998年全国大学生数学建模竞赛的A题问题如下：市场上有n种资产（股票、债券、）Si（i=1，n）供投资者选择，某公司有数额为M的一笔相当大的资金可用作一个时间的投资。公司财务分析人员对这n种资产进行了评估，估算出在这一时期内购买Si有平均收益率为ri，并预测出购买Si的风险损失率为qi。考虑到投资越分散总的风险越小，公司确定，当用这笔资金购买若干种资产时，总体风险可用所投资的Si中最大的一个风险来度量。购买Si要付交易费，费率为pi，并且当购买额不超过给定值ui时，交易费按购买ui计算（不买当然无须付费）。另外，假定同期银行存款利率是r0，且既无交易费又无风险（r0=5%）。（1）已知n=4时的相关数据如下：试给该公司设计一种投资组合方案即用给定的资金M，有选择地购买若干种资产或存银行生息，使净收益尽可能大，而总体风险尽可能小。 （2）试就一般情况对以上问题进行讨论，利用以下数据进行计算。 在AD边等距地设置7个波源Ri（i=1，7），BC边对等地安放7个接收器S j（j=1，7），记录由Ri发出的弹性波到达Sj的时间ij秒），见表2-3。 ij S1 S2 S3 S4 S5 S6 S7 R1 0.0645 0.0602 0.0813 0.3516 0.3867 0.4314 0.5721 R2 0.0753 0.0700 0.2852 0.4341 0.3491 0.4800 0.4980 R3 0.3456 0.3205 0.0974 0.4093 0.4240 0.4540 0.3112 R4 0.3655 0.3289 0.4247 0.1007 0.3249 0.2134 0.1017 R5 0.3165 0.2409 0.3214 0.3256 0.0904 0.1874 0.2130 R6 0.2749 0.3891 0.5895 0.3016 0.2058 0.0841 0.0706 R7 0.4434 0.4919 0.3904 0.0786 0.0709 0.0914 0.0583 2、0-1型变量在数学建模中的应用型变量在数学建模中的应用1、空洞探测问题、空洞探测问题1.1 问题的提出问题的提出 这是2000年全国大学生数学建模竞赛的D题。 山体、隧洞、坝体等的某些内部结构可用弹性波测量来确定。一个简化问题可描述为，一块均匀介质构成的矩形平板内有一些充满空气的空洞，在平板的两个邻边分别等距地设置若干波源，在它们的对边对等地安放同样多的接收器，记录弹性波由每个波源到达对边上每个接收器的时间，根据弹性波在介质中和空气中不同的传播速度，来确定板内空洞的位置。现考察如下的具体问题：一块240（米）240（米）的平板 （如图21）在AB边等距地设置7个波源Pi（i=1，7），CD边对等地安放7个接收器Q j（j=1，7），记录由Pi发出的弹性波到达Qj的时间tij（秒），见表2-2； tij Q1 Q2 Q3 Q4 Q5 Q6 Q7 P1 0.0611 0.0895 0.1996 0.2032 0.4181 0.4923 0.5646 P2 0.0989 0.0592 0.4413 0.4318 0.4770 0.5242 0.3805 P3 0.3052 0.4131 0.0598 0.4153 0.4156 0.3563 0.1919 P4 0.3221 0.4453 0.4040 0.0738 0.1789 0.0740 0.2122 P5 0.3490 0.4529 0.2263 0.1917 0.0839 0.1768 0.1810 P6 0.3807 0.3177 0.2364 0.3064 0.2217 0.0939 0.1031 P7 0.4311 0.3397 0.3566 0.1954 0.0760 0.0688 0.1042 在AD边等距地设置7个波源Ri（i=1，7），BC边对等地安放7个接收器S j（j=1，7），记录由Ri发出的弹性波到达Sj的时间ij秒），见表2-3。 ij S1 S2 S3 S4 S5 S6 S7 R1 0.0645 0.0602 0.0813 0.3516 0.3867 0.4314 0.5721 R2 0.0753 0.0700 0.2852 0.4341 0.3491 0.4800 0.4980 R3 0.3456 0.3205 0.0974 0.4093 0.4240 0.4540 0.3112 R4 0.3655 0.3289 0.4247 0.1007 0.3249 0.2134 0.1017 R5 0.3165 0.2409 0.3214 0.3256 0.0904 0.1874 0.2130 R6 0.2749 0.3891 0.5895 0.3016 0.2058 0.0841 0.0706 R7 0.4434 0.4919 0.3904 0.0786 0.0709 0.0914 0.0583 已知弹性波在介质和空气中的传播速度分别为2880（米/秒）和320（米/秒），且弹性波沿板边缘的传播速度与在介质中的传播速度相同。1）确定该平板内空洞的位置。2）只根据由Pi发出的弹性波到达Qj的时间tij（i，j=1，7），能确定空洞的位置吗？讨论在同样能够确定空洞位置的前提下，减少波源和接受器的方法。1.2 模型的假设及符号说明（1）模型的假设波源和接收器的设置仅按图中的设置来考虑，不考虑其它情况。在图中的一个小方格要么全是空气，要么全。 是介质。（2）符号说明符号说明n+1：x轴方向等距地放置波源Pi（i=17）的个 数，相应地n是x轴方向的方格数（已知 的n=6）；m+1：y轴方向等距地放置波源Ri（i=17）的个 数，相应地m是y轴方向的方格数（已知的 m=6）；d1： x轴方向方格的边长（已知的，d1=40米 ）；d2： y轴方向方格的边长（已知的，d2=40米 ）；tij（i=06；j=06）：由波源Pi发生的弹性波到达 接收器Qj的时间；pij（i=06；j=06）：波线PiQj经历的介质的 总长度；qij（i=06；j=06）：波线PiQj经历的空洞的 总长度；v1：弹性波在介质中的传播速度（已知的 v1=2880米/秒）；v2：弹性波在空气中的传播速度（已知的 v2=320米/秒）； 1.3 问题的分析及数学模型问题的分析及数学模型 （1）问题的分析这是一个反问题，不妨先看它正问题。已知空洞的位置及弹性波在介质和空气中的传播速度，求由波源Pi发出的弹性波到达接收器Qj的时间tij，为了求出tij，可先算出波线PiQj经历的每个小方格的长度，根据该方格是介质还是空洞，分别除以波在介质和空气中的传播速度，对经历的所有方格求和，即得tij。上述的正问题等价于：算出波线PiQj经历的各个小方格的长度，根据该方格是介质还是空洞，分别乘以0和1，对经历的所有方格求和，即可得到波线PiQj经历的空洞的长度。反问题可以这样求解：先由波源Pi发出的弹性波到达接收器Qj的时间tij及弹性波在介质和空气中的传播速度，算出波线PiQj经历的空洞的总长度，再由PiQj经历的每个方格的长度，即可解出每个方格是介质还是空洞（即是0还是1）。（2）数学建模）数学建模建立坐标系如图2-2，X、Y轴上分别等距地放置n+1，m+1个波源（本题n=m=6），将平板划分为 个方格，则1）波线PiQj经历的介质的总长度pij和空洞的总长度qij满足下面两个式子：)36(nmijijijtvqvp21xxyo126781231 3236345132）波线PiQj的方程为 (3）直线方程(3)与 的交点可算出,从而求得PiQj经历的每个方格长度。3）将 个方格如图示顺序排列，第i个方格的特征量记为 (4） 全体方格的特征向量记为 。), 1 , 0,()(njimiyjimx(0 ),(0 )x k kn y k km)36(nm个方格为空气第个方格为介质第iixi10Tmnxxx),(12221)(mddjiqpijij4）由AB至CD的波线PiQj共 条，其中n+1条PiQj是无用的，去掉无用波线后n(n+1）条波线PiQj的顺序按 ， ； , ； 排列，每条波线经历的各个方格的长度排成一行向量，于是得到 条波线经历mn个方格的长度矩阵 5）将1）得到的波线PiQj经历空洞的总长度qij，按照4）中波线的顺序构成空洞长度（列）向量 ，则应有 （5）2) 1( n0inj11i2 jn) 1(0,njni) 1( nn) 1 () 1(mnnnA)1 ()1( nnb)1()1(bxA0j 6）完全类似地处理波线RkSl，得到m（m+1）条波线经历mn个方格的长度矩阵 和空洞长度（列）向量 ，构造 , 则)2()1(mnmmA)2()1(mmb)2()1()1()1(AAAmnmmnn)2()1(bbbbAx 当 Rank(A)=mn 时可求出该方程最小二乘意义下的解bAAAxTT1)(7）取适当的 （8） 取 的小方格为介质， 小方格为空洞，其余的（如果有的话）无法确定。 2.3.4 模型求解模型求解% 空洞探测ch23.m% 文件名：ch23.m% 计算 PiQj 在每个格子中的线段长度。clearn=1;1110iiixxx令, 00ix1ixfor i=0:6 for j=0:6 ai=i;aj=j; if i=j aj=j+1; if aj=7 break; end endChudian=ai,0;Modian=aj,6;Dianwei=ai,0;pp=; for m=1:6 x=Chudian(1)+(Modian(1)-Chudian(1)*(m-Chudian(2)/. (Modian(2)-Chudian(2); Dianwei=Dianwei;x,m;end if ajai p=aj;aj=ai;ai=p;end if abs(aj-ai)=2 for k=(aj+1):(ai-1) y=Chudian(2)+(Modian(2)-Chudian(2)*(k-Chudian(1)/. (Modian(1)-Chudian(1); pp=pp;k,y;end end Jiaodian=Dianwei;pp;hang,lie=size(Jiaodian); aa=Jiaodian(:,1);bb=unique(aa);hang1,lie1=size(bb); Zuobiao=; for ii=1:hang1 for jj=1:hang if Jiaodian(jj,1)=bb(ii) Zuobiao(ii,:)=Jiaodian(jj,:);break;end;end;end; for pp=1:hang1-1 d=sqrt(sum(Zuobiao(pp,:)-Zuobiao(pp+1,:).2); if j=2 ifdistance_PQ(n,:)=distance_PQ(n-1,:); distance_PQ(n,:)=; else n=n+1; end else n=n+1; end endend%计算矩阵 A1 的秩A1=distance_PQ;disp(矩阵 A1 的秩)rank(A1)%计算RiSj在每个格子中的线段长度。n=1;for i=0:6 for j=0:6 ai=i;aj=j; if i=j aj=j+1; if aj=7break; end endChudian=0,ai;Modian=6,aj;Dianwei=0,ai;pp=; for m=1:6 y=Chudian(2)+(Modian(2)-Chudian(2)*(m-Chudian(1)/. (Modian(1)-Chudian(1); Dianwei=Dianwei;m,y;end if ajai p=aj;aj=ai;ai=p;end if abs(aj-ai)=2 for k=(aj+1):(ai-1) x=Chudian(1)+(Modian(1)-Chudian(1)*(k-Chudian(2)/. (Modian(2)-Chudian(2); pp=pp;x,k;end endJiaodian=Dianwei;pp;hang,lie=size(Jiaodian); aa=Jiaodian(:,2);bb=unique(aa);hang1,lie1=size(bb);Zuobiao=; for ii=1:hang1 for jj=1:hang if Jiaodian(jj,2)=bb(ii) Zuobiao(ii,:)=Jiaodian(jj,:);break;end;end;end; for pp=1:hang1-1 d=sqrt(sum(Zuobiao(pp,:)-Zuobiao(pp+1,:).2); if j=2 ifdistance_RS(n,:)=distance_RS(n-1,:); distance_RS(n,:)=; else n=n+1; end else n=n+1; end endend%对矩阵进行合并操作，计算系数矩阵 AA=distance_PQ;distance_RS;disp(矩阵 A 的秩)rank(A)%首先输入时间矩阵,并给出 bb1=0.0895 0.1996 0.2032 0.4181 0.4923 0.5646; 0.0989 0.4413 0.4318 0.4770 0.5242 0.3805; 0.3052 0.4132 0.4153 0.4156 0.3563 0.1919; 0.3221 0.4453 0.4040 0.1789 0.0740 0.2122; 0.3490 0.4529 0.2263 0.1917 0.1768 0.1810; 0.3807 0.3177 0.2364 0.3064 0.2217 0.1031; 0.4311 0.3397 0.3566 0.1954 0.0760 0.0688; b2=0.0602 0.0813 0.3516 0.3867 0.4314 0.5721; 0.0753 0.2852 0.4341 0.3491 0.4800 0.4980; 0.3456 0.3205 0.4093 0.4240 0.4540 0.3112; 0.3655 0.3289 0.4247 0.3249 0.2134 0.1017; 0.3165 0.2409 0.3214 0.3256 0.1874 0.2130; 0.2749 0.3891 0.5895 0.3016 0.2058 0.0706; 0.4434 0.4919 0.3904 0.0786 0.0709 0.0914;lie1=;lie2=;for i=1:7 m=b1(i,:);n=b2(i,:); lie1=lie1;m;lie2=lie2;n;endb=lie1;lie2;%计算方程 Ax=b 最小二乘意义下解x=Ab;tushi=;for i=0:5 pp=x(i*6+1:(i+1)*6);tushi=tushi;pp;endx=tushi; disp(输出结果 X )x%进行空洞判断a=0.11;b=-0.11;for i=1:6 for j=1:6 if tushi(i,j)b & tushi(i,j)a xx(i,j)=0; else xx(i,j)=1; end endendxx执行后输出矩阵 A1 的秩ans = 29矩阵 A 的秩ans = 36输出结果 X x = 0.0098 0.0180 0.0058 0.0070 0.0181 0.0137 0.0145 0.1233 0.1325 0.0140 0.0137 0.0125 0.0203 0.1196 0.1204 0.0217 0.1155 0.0174 0.0255 0.0128 0.1273 0.1246 0.0173 0.0074 0.0127 0.1271 0.0145 0.0134 0.0069 0.0134 0.0054 0.0159 0.0134 0.0075 0.0185 0.0173取=0.1有 xx = 0 0 0 0 0 0 0 1 1 0 0 0 0 1 1 0 1 0 0 0 1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0对应于平板上的空洞位置如图23（注意:相对于x的排列，上下正相反）。ABCD2.3.5 问题的进一步讨论1）只根据由Pi发出的弹性波到达Qj的时间tij，可得长度矩阵 ，若(5)式有解则能确定空洞的位置，但 ，(5)式无解，故不能确定空洞的位置。2）在同样能够确定空洞位置的前提下，减少波源和接受器的方法：设在BC边仍安放7个接收器Sl，而AD边设置波源Rk的数量一个一个地增加（由下而上），发现设置两个波源时有 ， 虽然（7）式有解，但由（7）式的解无法选取得到 。直到设置5个波源时由（7）式才得到满意的解。)1(3642A29)()1 (ARank36)(ARankxx= -0.0324 0.0485 -0.0889 -0.0959 0.0231 0.0661 0.0181 0.9766 1.0591 0.0345 0.0201 -.0688 0.0563 0.9402 0.9846 0.1003 0.9100 -0.0468 0.0588 0.0430 1.0429 0.9694 0.0085 -0.0198 0.1416 0.9766 -0.1041 -0.0225 -0.0726 -.0070 -0.2188 0.0021 0.0664 -0.0444 0.0680 0.0376取=0.22，即得到与前面同样的结果。 另外，也可试验在AB和AD边的各7个点上，选取不同位置设置波源，并类似地考虑接收器，找出减少波源（和接收器）的方案。 2.1.3 MATLAB求解对于方程组 ，若A是满秩方阵，则MATLAB求解命令为： MATLAB中常用的矩阵运算命令 1、求解方程组 若 A 是满秩的方阵，求解命令为：它等价于2、求解方程组 若 A 不是满秩的方阵，求解命令为：rrefC,其中C=A b3 、矩阵A的秩 rank（A）4、V，D=eig(A) 输出的D是矩阵A的特征值构成的对角矩阵，输出的V是矩阵A的特征向量构成的特征向量矩阵； bAx bAxbAx1bAx