运筹学实验报告.doc

精品文档 运筹与优化课内实验 检测 实验一：线性规划问题与对偶问题的建模与求解一 线性规划：满足以下三个条件的称之为线性规划问题：（1）决策变量的取值是连续的，既可以为整数，也可以喂分数、小数或实数。（2）目标函数是决策变量的线性函数。（3）结束条件是含决策变量的线性等式或者不等式。二线性规划模型的形式：2.1一般形式 (2.1)矩阵形式 (2.2)其中为决策向量，为目标函数的系数向量，为常数向量，为系数矩阵。2.2标准形式所谓线性规划问题的标准形式，是指目标函数要求所有约束条件都是等式约束，且所有决策定量都是非负的，即 三原问题与对偶问题的表达形式和关系 在线性规划的对偶理论中，把如下线性规划形式称为原问题的标准形式而把如下线性规划形式称为对偶问题的标准形式若用矩阵形式表示，则原问题和对偶问题分别可写成如下形式：原问题 对偶问题 原问题与对偶问题的关系见表原问题（对偶问题）对偶问题（原问题）minmax目标函数系数右边常数约束条件系数矩阵右边常数目标函数系数系数矩阵的转置第个约束条件为“”型第个约束条件为“=”型第个约束条件为“”型第个变量0第个变量无限制第个变量0第个变量0第个变量无限制第个变量0第个约束条件为“”型第个约束条件为“=”型第个约束条件为“”型四实际问题与lingo求解计算例4.1.1（原问题）：设某种植物每天至少需要2g水，4g矿物质，5g维生素。已知三种肥料A，B，C；其中A种肥料含有1g水，3g维生素；B种肥料含有3g水，2g矿物质，1g维生素；C种肥料含有1g水，2g矿物质，2g维生素。其中A种肥料6元每克，B种肥料4元每克，C种肥料7元每克，要求确定既要满足植物生长的营养需求，又要费用最省的选用肥料的方案。该问题的模型为:min z=6*x1+4*x2+7*x3; s.t.用lingo求解代码如下：model:min=6*x1+4*x2+7*x3;x1+3*x3=2;3*x1+2*x2+x3=4;x1+2*x2+2*x3=5;x1=0;x2=0;x3=0;end求解得：Global optimal solution found at iteration: 0 Objective value: 12.00000 Variable Value Reduced Cost X1 0.000000 3.000000 X2 1.833333 0.000000 X3 0.6666667 0.000000 Row Slack or Surplus Dual Price 1 12.00000 -1.000000 2 0.000000 -1.000000 3 0.3333333 0.000000 4 0.000000 -2.000000 5 0.000000 0.000000 6 1.833333 0.000000 7 0.6666667 0.000000通过上面结果可以知道，最佳费用为12元，此时，需要A,B,C三种肥料为0,1.833,0.667克。例4.1.2（对偶形式）：设某公司生产A，B，C三种肥料，已知生产A种肥料赢利2元，生产B种肥料赢利4元，生产C种肥料赢利5元。而生产A需要1小时设备1，3小时设备2，1小时设备3；生产B需要 2小时设备2，2小时设备3；需要3小时设备1，1小时设备2，2小时设备3；而设备1每天可用6小时，设备2每天可用4小时，设备3每天可用7小时，求怎样安排生产才能使公司的赢利最大。该问题的模型为： max z=2*y1+4*y2+7*y3;s.t.用lingo求解model:max=2*y1+4*y2+5*y3;y1+3*y2+y3=6;2*y2+2*y3=4;3*y1+y2+2*y3=0;y2=0;y3=0;end求解得： Global optimal solution found at iteration: 0 Objective value: 12.00000 Variable Value Reduced Cost Y1 1.000000 0.000000 Y2 0.000000 0.3333333 Y3 2.000000 0.000000 Row Slack or Surplus Dual Price 1 12.00000 1.000000 2 3.000000 0.000000 3 0.000000 1.833333 4 0.000000 0.6666667 5 1.000000 0.000000 6 0.000000 0.000000 7 2.000000 0.000000例4.2.1（原问题）：某企业利用3种原料B1，B2，B3生产2种产品A1，A2。3种原料的月供应量和生产1吨产品A1，A2均可在市场销售，企业应如何安排月生产计划，使总收益最大？单位产品消耗 原料A1A2原料月供应量（t）B1B2B3123132150240300单位产品价格（万元/t）2.41.8 设计划生产产品Ai的数量为xi（t/月），i=1，2，则所讨论的原问题的数学模型为：max z=2.4*x1+1.8*x2;s.t. 用lingo求解model:max=2.4*x1+1.8*x2;x1+x2=150;2*x1+3*x2=240;3*x1+2*x2=0;x2=0;end求解得：Global optimal solution found at iteration: 0 Objective value: 244.8000 Variable Value Reduced Cost X1 84.00000 0.000000 X2 24.00000 0.000000 Row Slack or Surplus Dual Price 1 244.8000 1.000000 2 42.00000 0.000000 3 0.000000 0.1200000 4 0.000000 0.7200000 5 84.00000 0.000000 6 24.00000 0.0000004.2.2（对偶问题）：设原料B1，B2，B3的价格为y1,y2,y3(万元/t)，显然，应有yi=0,I=1,2,3.由原问题的条件，生产1t产品A1需消耗1t原料B1，2t原料B2,3t原料B3，可获得收益为2.4万元。因此，若将生产1t产品A1的这些原料卖出所得的收益为y1+2*y2+3*y3 (万元)它必须不少于生产1t产品A1所得的收益，对于该企业才是合算的。所以应有y1+2*y2+3*y3=2.对于产品A2，可类似得到y1+2*y2+3*y3=1.8。同时，若买方（公司）欲购买该工厂的全部原料，则应付出150*y1+240*y2+300*y3万元（这也是该工厂卖出全部原料的总收益）。从买方角度应使总支出金可能地少。因此，可得到线性规划问题。该问题的模型为:min f=150*y1+240*y2+300*y3s.t. 用lingo求解model:min=150*y1+240*y2+300*y3；y1+2*y2+3*y3=2.4y1+3*y2+2*y3=1.8;y1=0;y2=0;y3=0;求解得： Global optimal solution found at iteration: 4 Objective value: 244.8000 Variable Value Reduced Cost Y1 0.000000 42.00000 Y2 0.1200000 0.000000 Y3 0.7200000 0.000000 Row Slack or Surplus Dual Price 1 244.8000 -1.000000 2 0.000000 -84.00000 3 0.000000 -24.00000 4 0.000000 0.000000 5 0.1200000 0.000000 6 0.7200000 0.000000五实验分析1.针对上面的线性规划问题，每个问题都相应的给出了具体的对偶问题，方便的比较。2.利用lingo软件求解，速度快，精度高。实验二：基于匈牙利解法的指派问题建模与求解一、指派问题与匈牙利解法指派问题又称分配问题，其用途非常广泛，比如某公司指派n个人去做n件事，各人做不同的一件事，如何安排人员使得总费用最少？若考虑每个职工对工作的效率（如熟练程度等），怎样安排会使总效率达到最大？这些都是一个企业经营管理者必须考虑的问题，所以该问题有重要的应用价值.虽然指派问题可以用0-1规划问题来解，设X(I,J)是0-1变量, 用X(I,J)=1表示第I个人做第J件事, X(I,J)=0表示第I个人不做第J件事. 设非负矩阵C（I，J）表示第I个人做第J件事的费用， 则问题可以写成LINGO程序SETS:PERSON/1.N/;WORK/1.N/;WEIGHT(PERSON, WORK): C, X ;ENDSETSDATA:W=ENDDATAMIN= SUM(WEIGHT: C*X);FOR(PERSON(I): SUM(WORK(J):X(I,J)=1);FOR(WORK(J): SUM(PERSONM(I):X(I,J)=1);FOR(WEIGHT: BIN(X);其中2*N个约束条件是线性相关的, 可以去掉任意一个而得到线性无关条件.但是由于有N2个0-1变量, 当N很大时，用完全枚举法解题几乎是不可能的. 而已有的0-1规划都是用隐枚举法做的，计算量较大. 对于指派问题这种特殊的0-1规划，有一个有效的方法匈牙利算法，是1955年W. W. Kuhn利用匈牙利数学家D.Knig的二部图G的最大匹配的大小等于G的最小顶点覆盖的大小的定理提出的一种算法，这种算法是多项式算法，计算量为O(N3).匈牙利算法的基本原理是基于以下两个定理. 定理1 设C=(Cij)nn是指派问题的效益矩阵，若将C中的任一行（或任一列）减去该行（或该列）中的最小元素，得到新的效率矩阵C，则C对应的新的指派问题与原指派问题有相同的最优解. 定理2 效率矩阵C中独立的0元素的最多个数等于覆盖所有0元素的最少直线数. 当独立零元素的个数等于矩阵的阶数时就得到最优解. 二模型扩展1、目标函数最大化的指派问题模型对于最大化的指派问题 不能将目标函数改为来实现求解，因为匈牙利算法要每一个系数都为非负. 因此，对于最大化的指派问题，只能在效率矩阵C=(cij)nn找出最大的数M，然后进行矩阵变换，即可以讨论矩阵B=M C, 其中M=max(max (C)*ones( size(C)，将原指派问题模型转化为同解指派问题模型： 再用匈牙利算法求解即可.2、效率矩阵不是方阵的指派问题模型若人多事少, 可以添上一些虚拟的事, 这些事被各人做的费用可取为零.若人少事多, 可以添上一些虚拟的人, 他们做各事的费用可取为零.若一个人可同时做几件事, 可以将这人化为相同的几个人, 费用不变.若规定某人不能做某事, 求最小时则可令这人做这事的费用为充分大的数(在MATLAB中可取为inf ) ；求最大时则可令这人做这事的费用为0，然后按目标函数最大化的指派问题数学模型的求解方式求解. 三匈牙利解法的缺陷与改进3.1、算法的不足之处指派问题就是给定一个n阶非负矩阵A=A(I,J), 求一个0-1 矩阵X(I,J), 使得众所周知, 求这个问题的一个方法是W.W.Kuhn在1955年利用匈牙利数学家D.Knig的关于矩阵中独立0元素的定理(即矩阵中独立零元素的最大个数等于覆盖全部0元素的直线的最少条数) 提出了解指派问题的一种算法. 步骤一: 将系数矩阵A的各行元素分别减去本行中的最小元素, 再将所得矩阵(仍记为A)的各列元素分别减去本列中的最小元素. 这样所得矩阵A的各行和各列中都有0元素.步骤二: 求A的最大个数的独立0元素，如个数等于n， 则完成. 不然进行下一步.步骤三: 求覆盖矩阵A所有0元素的最少数目的直线, 求A中所有未被直线覆盖的元素中的最小数m， 把未被直线覆盖的行中各元素都减去m，再把A的被直线覆盖的各列中各元素都加上m. 转步骤二.这是一种正确的多项式算法，但是我们发现现有的文献中具体实现步骤二的算法是有缺陷的，叙述如下：现有的求最大个数的独立0元素的算法是:步骤一: 找只有一个没有标记的0元素的行, 如无转步骤三步骤二: 将该0元素标记为独立0元素, 并把该独立0元素所在的列的其他0元素标记为非独立0元素, 转步骤一步骤三: 找只有一个没有标记的0元素的列, 如无转步骤五步骤四: 将该0元素标记为独立0元素, 并把该独立0元素所在的行的其他0元素标记为非独立0元素, 转步骤三步骤五: 找没有标记的0元素, 如无, 完成最大个数的独立0元素的寻求. 步骤六: 未标记的0元素在每行及每列的个数都大于1，按某种规则取某个0元素标记为独立0元素, 并把该独立0元素所在的行和列的其他0元素标记为非独立0元素, 转步骤一.对于已求出最大个数的独立0元素但个数小于n的矩阵A, 一般文献上介绍的求覆盖全部零元素的方法是：步骤一: 对所有无独立0元素的行打;步骤二: 若无新打的行, 转步骤六;步骤三: 对新打的行中的非独立0元素所在的列打;步骤四: 若无新打的列, 转步骤六;步骤五: 对新打 的列中的独立0元素所在的行打. 转步骤二;步骤六: 用直线覆盖没有打的行与打的列.3.2、算法漏洞分析但是我们发现, 在求最大个数的独立0元素的方法中的步骤六有漏洞. 就是说在当行和列中未标记的0元素的个数都大于1时在独立0元素的确定上可能会把不应当是独立0元素的0元素错误地确定为独立0元素，从而得不到最大个数的独立0元素. 一般文献上介绍的求最大个数的独立0元素的方法中，当矩阵中各行各列的无标记的0元素的个数都大于1时，求独立0元素的规则有两种:规则1: 可任选其中一个0元素为独立0元素（本文用上标加撇来标记：0），同时标记同行和同列中其它0元素为非独立0元素（本文用加一横线来标记：）. 规则2: 若各行中无标记的0元素至少有两个，. 从剩有0元素最少的行开始，比较这行各0元素所在列中0元素的数目，选择0元素少的那列的这个0元素为独立0元素. 然后标记同行同列的其它0元素为非独立0元素.现有的算法中的缺点在于在确定独立零元素时不一定能得到最大个数的独立零元素，从而算法会失败. 我们找到以下一个有6个独立0元素的6阶矩阵A, 其中每行每列都有二个或2个以上的0元素，其中用0表示是独立0元素，表示划线的0元素，1表示任意正数):但是按照取第一行第一列的0元素为独立0元素时只能得到以下所示5个独立0元素，小于最大独立0元素的个数6，而需要覆盖全部零元素的直线的条数是6条竖线，导致算法失败：从以上例子可见文献中叙述的独立0元素的求法不一定能得到最大数目的独立零元素. 我们可以用两种方法构造任意大于6阶的类似的例子，如以下10阶的两个例子分别应该有10个独立0元素，但按不当地选取独立0元素时只能得到9个独立0元素.通过把这些矩阵作为对角块，其余元素为正数的矩阵可以构造更多的这种例子.3.3、算法的更正由于问题出在求最大个数的独立零元素上，我们在算法中增加了利用求M-增广路的方法3来保证能得到最大个数的独立零元素.我们提出的更正的完整的算法如下:一、化约矩阵将系数矩阵A的各行元素分别减去本行中的最小元素, 再将所得矩阵(仍记为A)的各列元素分别减去本列中的最小元素. 这样所得矩阵A的各行和各列中都有0元素.二、标记零元素 找只有一个没有标记的0元素的行, 将该0元素标记为0(独立0元素), 并把该0所在的列的其他0元素标记为(划线零元素), 转1. 如无转2. 找只有一个没有标记的0元素的列, 将该0元素标记为0, 并把该0所在的行的其他0元素标记为, 转2；如无转3. 若未标记的0元素在其所在的行及列的个数都大于1，取未标记的0元素最少的行或列中取一个所在列或行上0元素个数较少的零元素标记为0, 并把该0所在的行和列的其他0元素标记为, 转1；如无,转三.三、求最小点覆盖(即最小直线覆盖) 若0的个数等于n，完成；否则，转2. 对所有没有0的行标记为S, 对S行中的所在列标记为T, 将T列上0所在的行标记为TS，将S行T列上标记为. 将S行和T列分别标记为S和T，将TS改为S. 如S行中的某行上的位于未标志的列上，标记该列为T，如所有的T列上都有0, 将0所在的行标为TS. 将S行T列上的标为为, 转2. 不然，若S行中的某行上的位于未标志某个列上，而该列上没有0，则存在M-增广路，在增广路上交换和0. 取消行、列的标记，转三；若S行上不存在位于未标志的列上的, 将S改为S，转四.四、进一步化约矩阵用直线覆盖标记为T的列和没有标记为S的行. 即最小点覆盖是(X-S)T, 求矩阵中所有未被直线覆盖的元素中的最小数m， 把未被直线覆盖的行中各元素都减去m，再把被直线覆盖的各列中各元素都加上m. 去掉所有的标记，转二. 四实际问题分析与求解4.1（基于匈牙利解法的问题求解）：某出版社有A、B、C、D、E五名专业翻译，他们翻译五种外文的速度（印刷字符/小时）如表所示.若规定一名翻译只能翻译一种外文，问：若C翻译不懂法语、D翻译不懂英语（如表所示）如何安排效率最高？若根据原文作者的要求，A翻译必须翻译德语，E翻译必须翻译日语，如何安排效率最高？语种翻译英法俄德日A1000700800800600B80060010009001100C900/700900900D/1000700700800E700700700900900解析：此问题是一个有条件的指派问题的求解，C翻译不懂法语、D翻译不懂英语，故设翻译C翻译法语、翻译D翻译英语的效率为零.于是，该有条件的指派问题的效率矩阵C为：目标函数为：.将其化为求极小化问题：bij=maxcij-cij，于是所以，效率最高的最优安排方案是A翻译英语、B翻译俄语、C翻译德语、D翻译法语、E翻译日语.若根据原文作者的要求，A翻译必须翻译俄语、E翻译必须翻译日语，可视为此时只有B、C、D三名翻译需要安排任务，使之总的效率最高。语种翻译英法俄B8006001000C900/700D/1000700于是，该指派问题的效率矩阵C为：目标函数为：.将其化为求极小化问题：bij=maxcij-cij，于是所以，在满足A翻译必须翻译德语、E翻译必须翻译日语的前提下，B翻译俄语、C翻译英语、D翻译法语，这样分配时，效率最高.4.2（标准指派问题软件求解）：分配甲、乙、丙、丁、戊完成A，B，C，D，E五项任务，每人完成一项，每项任务只能由一个人完成，五个人完成任务所需时间如下表，做出合理安排，使总时间最少.人ABCDE甲8610912乙9127119丙74358丁958118戊467511（1）Lingo程序求解代码如下：model:sets:a/1.5/;time(a,a):t,n;endsetsdata:t=8 6 10 9 12 9 12 7 11 9 7 4 3 5 8 9 5 8 11 8 4 6 7 5 11;enddatamin=sum(time:t*n);for(a(i):sum(a(j):n(i,j)=1);for(a(j):sum(a(i):n(i,j)=1);end结果：Global optimal solution found at iteration: 6 Objective value: 30.00000 Variable Value Reduced Cost T( 1, 1) 8.000000 0.000000 T( 1, 2) 6.000000 0.000000 T( 1, 3) 10.00000 0.000000 T( 1, 4) 9.000000 0.000000 T( 1, 5) 12.00000 0.000000 T( 2, 1) 9.000000 0.000000 T( 2, 2) 12.00000 0.000000 T( 2, 3) 7.000000 0.000000 T( 2, 4) 11.00000 0.000000 T( 2, 5) 9.000000 0.000000 T( 3, 1) 7.000000 0.000000 T( 3, 2) 4.000000 0.000000 T( 3, 3) 3.000000 0.000000 T( 3, 4) 5.000000 0.000000 T( 3, 5) 8.000000 0.000000 T( 4, 1) 9.000000 0.000000 T( 4, 2) 5.000000 0.000000 T( 4, 3) 8.000000 0.000000 T( 4, 4) 11.00000 0.000000 T( 4, 5) 8.000000 0.000000 T( 5, 1) 4.000000 0.000000 T( 5, 2) 6.000000 0.000000 T( 5, 3) 7.000000 0.000000 T( 5, 4) 5.000000 0.000000 T( 5, 5) 11.00000 0.000000 N( 1, 1) 0.000000 0.000000 N( 1, 2) 0.000000 0.000000 N( 1, 3) 0.000000 3.000000 N( 1, 4) 1.000000 0.000000 N( 1, 5) 0.000000 3.000000 N( 2, 1) 0.000000 1.000000 N( 2, 2) 0.000000 6.000000 N( 2, 3) 0.000000 0.000000 N( 2, 4) 0.000000 2.000000 N( 2, 5) 1.000000 0.000000 N( 3, 1) 0.000000 3.000000 N( 3, 2) 0.000000 2.000000 N( 3, 3) 1.000000 0.000000 N( 3, 4) 0.000000 0.000000 N( 3, 5) 0.000000 3.000000 N( 4, 1) 0.000000 2.000000 N( 4, 2) 1.000000 0.000000 N( 4, 3) 0.000000 2.000000 N( 4, 4) 0.000000 3.000000 N( 4, 5) 0.000000 0.000000 N( 5, 1) 1.000000 0.000000 N( 5, 2) 0.000000 4.000000 N( 5, 3) 0.000000 4.000000 N( 5, 4) 0.000000 0.000000 N( 5, 5) 0.000000 6.000000 Row Slack or Surplus Dual Price 1 30.00000 -1.000000 2 0.000000 -9.000000 3 0.000000 -9.000000 4 0.000000 -5.000000 5 0.000000 -8.000000 6 0.000000 -5.000000 7 0.000000 1.000000 8 0.000000 3.000000 9 0.000000 2.000000 10 0.000000 0.000000 11 0.000000 0.000000（2）MATLAB求解此模型的代码:C=8 6 10 9 12;9 12 7 11 9;7 4 3 5 8;9 5 8 11 8;4 6 7 5 11;n=size(C,1);C=C(:);A=;B=;Ae=zeros(2*n,n2);for i=1:n for j=(i-1)*n+1:n*i Ae(i,j)=1; end for k=i:n:n2 Ae(n+i,k)=1; endendBe=ones(2*n,1);Xm=zeros(n2,1);XM=ones(n2,1);x,z=linprog(C,A,B,Ae,Be,Xm,XM);x=reshape(x,n,n);disp(:);Assignment=round(x)disp(:);zOptimization terminated.最优解矩阵为:Assignment = 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0最优解为:z = 30.0000（3）基于改进匈牙利算法的matlab程序代码%主程序function f,D,G=assign(W)% 指派问题求最小费用的主程序，如求最大效益，改W=max(W(:)-W, f=max(W(:)*n-f即可% 输入- W: 指派问题的系数矩阵, W(I,J) 表示第I人做第j事的费用% 输出- D: 行向量, 第I人做第D(I)件事.% - G: 行向量，第J事让第G(J)人做. % - f: 最小费用null=min(W,2); f=sum(null);%null:各行最小元素; f:目标函数值的改变.n=numel(null); % n 权阵的阶数W=W-repmat(null,1,n);% 权阵各行减去本行最小元素null=min(W,1); f=f+sum(null);%null:权阵各列最小元素; f:目标函数值的改变.W=W-repmat(null,n,1);% 权阵各列减去本列最小元素，至此权阵各行各列至少有一个0while 1 D,G,E=assignz(W,n);% 调用标记独立0的程序assignz; S=find(D=0); % 无独立0的行标 if isempty(S) return; end; % 如有n个独立0, 完成指派工作 D,G,SP,TP=assignln(D,G,E,S,n); % 调用求最少划线程序assignln num=numel(SP)-numel(TP);% 最少划线数=n-num if num % 划线数小于n时, 变换权阵 null=setdiff(1:n,TP); % 此null是没有划线的列标 m=min(min(W(SP,null); % 此m为未被直线覆盖的元素中的最小数 W(SP,null)=W(SP,null)-m; % 未被直线覆盖的元素都减去m null=setdiff(1:n,SP); % 此add 是划过线的行标,标记独立0 W(null,TP)=W(null,TP)+m; % 余下的行和列加上 m f=f+m*num; % 目标函数值的改变. else return; % 划线数=n时完成end;end;%子程序function D,G,E=assignz(W,n) % 选独立0的子程序C=zeros(n);C(find(W=0)=1;E=C;%C:未标记的0;E:非独立的0D=zeros(1,n); G=zeros(1,n);%第I人做第D(I)事.第J事让第G(J)人做. u=sum(C,2); v=sum(C);%u(i)第i行中的未标记的0的个数,v(j)第j列中未标记的0的个数while any(u) row=find(u=1,1); % 行中只有一个未标记的0的第一个行标 while row col=find(C(row,:)=1,1); % 第row行0所在的列标 D(row)=col; G(col)=row; E(row,col)=0; u=u-C(:,col); v(col)=0; C(:,col)=0; % 删去第col列的0 row=find(u=1,1); % 为循环作预处理 end; % 这时行中不是只有一个未标记的0