前言MATLAB作为一款科学计算软件逐渐被广大科研人员所接受,以其强大的数据计算功能、图像的可视化界面及代码的可移植性受到了广大高校师生及科研人员的认可,借助MATLAB能够解决绝大部分的工程问题。对于从事数据分析和计算方面的工作者和学习者来说,MATLAB是一个很好的工具。MATLAB的创始人是Cleve Moler,他是美国工程院院士,MathWroks董事长和首席数学家。MATLAB可用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中,为需要进行有效数值计算的众多科学领域提供了一种全面的解决方案,并在很大程度上摆脱了传统非交互式程序设计语言如C、Fortran的编辑模式,代表了当今国际科学计算软件的先进水平。最优化理论和方法自古就有,最典型的实例就是黄金分割,最优化方法形成为科学方法则是在17世纪后。牛顿和莱布尼茨创立的微积分理论为最优化问题的解决提供了理论基础,而后产生的优化方法和简单库存模型等精确的解析方法可以称为古典最优化方法。虽然古典最优化方法出现得比较早,但是由于计算手段的限制,这些方法在解决实际问题时遇到了瓶颈,随着计算机的兴起和日益普及,原来制约优化设计的技术突飞猛进,得到了广泛的应用并创造了巨大的经济效益和社会价值。随着计算机应用技术的发展,各种可用于最优化的方法设计与实现的软件层出不穷,丰富了科研技术人员的开发研究手段,提高了解决实际问题的效率。而MATLAB以其强大的科学计算功能和覆盖面广、专业性强的工具箱发展成适合多学科、跨平台的大型实用科学计算软件,也为最优化计算问题的解决提供了有力的工具。本书以工程应用为基础,将最优化理论和方法与MATLAB相结合,帮助读者从理论和实践两个方面提高解决最优化的能力,让即便是数学基础不够深厚的读者,也同样能够利用MATLAB解决较难的最优化数学问题,为读者能够快速地进入这个领域、设计高效可行的最优化方案奠定一个扎实的基础。编写本书具有如下特点。1 版本新,函数新。MATLAB每年更新两次,神经网络工具箱也随之更新换代,许多旧的函数废弃不用,同时又有新的函数补充进来。本书基于MATLAB R2015b,介绍了新版本下的神经网络工具箱的使用方法。2 由浅入深,层次分明。本书的内容以最优化理论为主线,最优化方法与实际应用相结合的实例为基础,结合编者的多年教学实践经验,由浅入深地介绍各种最优化理论和方法在MATLAB中的实现方法。3 内容讲解不枯燥。本书结合相关理论和实践,由实践支撑理论,通过求解流程以及算法迭代过程的实现,让读者更容易理解并且掌握,书中许多实例是读者经常碰到的,读起来不枯燥。4 应用性强。在介绍最优化设计的章节中,每章都有应用MATLAB解决各领域中的实际最优化问题,最后以一个或几个实际应用的例子总结本章内容,帮助降低读者学习门槛、提高学习效率。全书共分为10章,主要内容包括:第1章介绍MATLAB R2015b,主要包括MATLAB的发展史、MATLAB R2015b的新功能、MATLAB的工作环境及控制流等内容。第2章介绍MATLAB计算基础,主要包括MATLAB的矩阵、绘图、编程技巧等内容。第3章介绍MATLAB数值计算,主要包括数据排序、符号运算、多项式运算、数据插值等内容。第4章对最优化计算进行概述,主要包括最优化的发展史、最优化的定义、优化工具箱等内容。第5章介绍线性规划,主要包括对线性规划的概述,线性规划的标准型、方法、实际应用等内容。第6章介绍非线性规划,主要包括非线性规划的概述、一维最优化方法、多维无约束非线性、非线性规划的实际应用等内容。第7章介绍整数规划,主要包括整数规划的概述、案例分析、求解、实际应用等内容。第8章介绍二次规划,主要包括等式约束二次规划及二次规划的MATLAB实现等内容。第9章介绍多目标规划及其他规划,主要包括多目标规划、最大最小化、半无限规划、动态规划等内容。第10章介绍了群智能算法,主要包括粒子群算法、遗传算法、模拟退火算法等内容。本书适用于计算应用、最优化领域和科学计算方向的研究人员使用,也可作为高校该类课程的本科生和研究生的教材,还可作为读者查询最优化数学问题求解方法的参考书。本书主要由李娅编写,此外参与编写的还有李炳辉、李丹、曾虹雁、邓俊辉、邓秀乾、邓耀隆、高泳崇、李嘉乐、李旭波、梁朗星、梁志成、刘超、刘泳、卢佳华、张棣华、张金林、钟东山、詹锦超、叶利辉、杨平和许兴杰。由于时间仓促,加之作者水平有限,错误和疏漏之处在所难免。在此,诚恳地期望得到各领域的专家和广大读者的批评指正。作者2016年5月
第5章 线性规划线性规划是运筹学中研究较早、发展较快、应用广泛、方法较成熟的一个重要分支,它是辅助人们进行科学管理的一种数学方法。是研究线性约束条件下线性目标函数的极值问题的数学理论和方法,英文缩写LP。它是运筹学的一个重要分支,广泛应用于军事作战、经济分析、经营管理和工程技术等方面。为合理地利用有限的人力、物力、财力等资源作出最优决策,提供科学的依据。5.1概述在经济管理、交通运输、工农业生产等经济活动中,提高经济效果是人们不可缺少的要求,而提高经济效果一般通过两种途径实现: 一是技术方面的改进,例如改善生产工艺,使用新设备和新型原材料; 二是生产组织与计划的改进,即合理安排人力、物力资源。线性规划研究的是: 在一定条件下,合理安排人力、物力等资源,使经济效果达到最好。一般的,求线性目标函数在线性约束条件下的最大值或最小值的问题,统称为线性约束条件的解,叫做可行解,由所有可行解组成的集合叫做可行域。决策变量、约束条件、目标函数是线性规划的三要素。5.1.1线性规划的发展史自从1947年乔治丹齐格(G.B.Dantzing)提出求解线性规划的单纯方法以来,线性规划在理论上趋向成熟,在应用上也日益深入,成为科学工程领域广泛应用的数学模型,其发展史如下。法国数学家J.B.J.傅里叶和C.瓦莱.普森分别于1832年和1911年独立地提出线性规划的想法,但未引起注意。1939年,苏联数学家Л.B.康托罗维奇在《生产组织与计划中的数学方法》一书中提出线性规划问题,也未引起重视。1947年,美国数学家G.B.Dantzing提出求解线性规划的单纯形法,为这门学科奠定了基础。1947年,美国数学家J.von诺伊曼提出对偶理论,开创了线性规划许多新的研究领域,扩大了它的应用范围和解题能力。1951年,美国经济学家T.C.库普曼斯把线性规划应用到经济领域,为此与康托罗维奇一起获得1975年的诺贝尔经济学奖。20世纪50年代后对线性规划进行了大量的理论研究,并涌现出一大批新的算法。例如,1954年C.莱姆基提出对偶单纯形法,1954年S.加斯和T.萨迪等人解决了线性规划的灵敏度分析和参数规划问题,1956年A.塔克提出互补松弛定理,1960年G.B.丹齐克和P.沃尔夫提出分解算法等。线性规划的研究成果还直接推动了其他数学规划问题(包括整数规划、随机规划和非线性规划的算法)的研究。由于数字电子计算机的发展,出现了许多线性规划软件,如MPSX、OPHEIE、UMPIRE等,可以很方便地求解几千个变量的线性规划问题。1979年,苏联数学家L.G.Khachian提出解线性规划问题的椭球算法,并证明它是多项式时间算法。1984年,美国贝尔电话实验室的印度数学家N.卡马卡提出解线性规划问题的新的多项式时间算法。用这种方法求解线性规划问题,在变量个数为5000时只要单纯形法所用时间的150。现已形成线性规划多项式算法理论。20世纪50年代后线性规划的应用范围不断扩大。5.1.2线性规划的一般问题下面探讨一个有关运输问题的线性规划实例,通过该实例将继续有关线性规划标准型的学习。【例51】(运输问题)设有两个建材厂C1和C2,每年沙石的产量分别为35万吨和55万吨,这些沙石需要供应到W1、W2和W3三个建筑工地,每个建筑工地对沙石的需求量分别为26万吨、38万吨和26万吨,各建材厂到建筑工地之间的运费(万元万吨)如表51所示,问题是: 应当怎样调运才能使得总运费最少?表51建材厂到建筑工地之间的运费工地运费万元万吨建材厂W1W2W3C110129C281113假设xij代表建材厂Ci运往建筑工地Wj的数量(万吨),则各建材厂和工地之间的运量可以用表52来表示。表52建材运送分配方案工地分配量万吨建材厂W1W2W3输出总量C1x11x12x1335C2x21x22x2355接收总量26382690可以得出总运费f=10x11 12x12 9x13 8x21 11x22 13x23。由表52可知,从建材厂C1运往各个工地的沙石的输出量应为建材厂C1的产量,即x11 x12 x13=35; 同样,从建材厂C1运往各个工地的沙石的输出量应为建材厂C2的产量,即x21 x22 x23=55。同时,各个工地的沙石需求量应当为从各个建材厂接收到的沙石的总量,因而有下面的约束条件:x11 x21=26x12 x22=38x13 x23=26注意到,建材厂要么给工地运送沙石,要么就不运送,因此运输量xij为一个非负值,即满足xij0。综合以上分析,整个运输问题可改成:min f=f=10x11 12x12 9x13 8x21 11x22 13x23x11 x12 x13=35x21 x22 x23=55x11 x21=26x12 x22=38x13 x23=26xij0(i=1,2; j=1,2,3)即求取满足约束条件的设计变量xij的值,使得总运费f最小。线性规划问题的具体情景虽然各有不同,但是在数学上均有以下共同特点:(1) 均可以用一组设计变量来表示一种实施方案。(2) 每个问题都有一定的约束条件,并可以用一组线性等式或者线性不等式表达。(3) 一般都有一个目标函数,该函数用于衡量方案的优劣,可以表达为设计变量的一个线性函数,目的一般为使得目标函数达到最大值或者最小值。5.2线性规划的标准型线性规划的目标函数可以是求最大值,也可以是求最小值,约束条件的不等号可以是小于号也可以是大于号。为了避免这种多样性形式带来的不便,一般采用统一的标准型进行描述。而如果遇到非标准型的线性规划问题,都可采用相应的方法将其改写成与其等价的线性规划标准型。5.2.1一般标准型根据线性规划问题的定义,线性规划问题即求取设计变量x=[x1,x2,,xn]T值,在线性约束条件下使得线性目标函数达到最大。由此可得,线性规划问题的一般标准形式如下:maxf=c1x1 c2x2 cnxns.t.a11x1 a12x2 a1nxn=b1a21x1 a22x2 a2nxn=b2am1x1 am2x2 amnxn=bmxi0i=1,2,,n其中,cii=1,2,,n,aiji=1,2,,m; j=1,2,,n,bii=1,2,,m均为给定的常数。5.2.2矩阵标准型利用向量或矩阵符号,线性规划问题的标准型还可以用矩阵形式表示:min f=cTxs.t.Axbx0通常A=aijmnRmn为约束矩阵,c=c1,c2,,cnTRn为目标函数系数矩阵; b=b1,b2,,bmTRm称为资源系数向量,x=x1,x2,,xnTRn称为决策向量。5.2.3向量标准型有时还将线性规划问题用向量的形式表示,此时线性规划的向量标准型为max f=cxs.t.[P1,P2,,Pn]x=bx0其中,Pj为矩阵A的第j列向量,例如:Pj=a1j,a2j,,amjT5.2.4非标准型根据实际应用问题建立的线性规划模型在形式上未必是标准型,对于不同类型的非标准型,可以采用相应的方法,通过以下方式将所建立的模型转化为线性规划的标准型。1. 目标函数为极小化设原有线性规划问题为极小化目标函数:min f=c1x1 c2x2 cnxn此时,可设f=-f,则极小化目标函数问题转化为极大化目标函数问题,即如下所示:max f=-c1x1 c2x2 cnxn2. 约束条件为不等式如果原有线性规划问题的约束条件为不等式,则可增加或减去一个非负变量,使约束条件变为等式,增加或减去的非负变量称为松弛变量。例如,约束为: ai1x1 ai2x2 ainxnbi可在左边增加一个非负变量xn 1,使其变为等式: ai1x1 ai2x2 ainxn xn 1=bi如果约束为: ai1x1 ai2x2 ainxnbi可在左边减去一个非负变量xn 1,使其变为等式: ai1x1 ai2x2 ainxn-xn 1=bi3. 无非负限制如果对某个变量xj的非负并无限制,可设两个非负变量xj和xj,令xj=xj-xj注意: 因为对原设计变量进行了代换,还需要将上式代入目标函数和其他约束条件进行相应的代换,这样即可满足线性规划标准型对变量非负的要求。【例52】将下列线性规划模型标准化,min f=x1-2x2 x3s.t.x1 x2 x35x1 x2-2x32-x1 2x2 3x3=6x1,x20将上述问题转化为等价的线性规划标准型。原问题的目标函数为求极小,即将目标函数两边乘以-1转化为求极大值,即求解目标为maxf=-x1 2x2-x3原问题约束条件中的前两个条件均为不等式,在第一个不等式的左边加上一个松弛变量x4,在第二个不等式的左边减去一个松弛变量x5,将两者转化为等式约束,即有x1 x2 x3 x4=5x1 x2-2x3-x5=2原问题对设计变量x2没有非负限制,因此,可引入非负变量x13和x23,令x3=x13-x23并将上式代入到目标函数和各约束条件中,最后整理可得与原问题等价的线性规划的标准型为max f=-x1 2x2-x3s.t.x1 x2 x3 x4=5x1 x2-2x3-x5=2-x1 2x2 3x13-3x23=6x1,x2,x13,x23,x4,x505.3单纯形法当线性规划的决策变量只有两个时,可以用图解法求解,但是一般线性规划问题的解法,常用的是单纯形法。对于标准形式的线性规划问题min fx=cTxs.t.Axbx0如果有有限最优值,则目标函数的最优值必在某一基本可行解处达到,因而只需在基本可行解中寻找最优解。这就有可能用穷举法来求得线性规划问题的最优解,但当变量很多时计算量很大,有时行不通。单纯形法(Simplex Method)的基本思想就是先找到一个基本可行解,检验是否为最优解或判断问题无解。否则,再转换到另一个使目标函数值减小的基可行解上,重复上述过程,直到求到问题的最优解或指出问题无解为止。设找到初始基本可行解x*,可行基为B,非基矩阵为N,即可写A=B,N。于是x*=B-1b0=b*0。相应的,目标函数值为f*=cx*=cB,cNb*0=cBb*,其中cB是c中与基变量xB对应的分量组成的m维行向量。再设任意可行解x=xBxN,由Ax=b得xB=B-1b-B-1NxN=b*-B-1NxN相应的目标函数值为f=cx=cBb*-cBB-1N-cNxN若记A=a1,a2,,an,则有f-f*-jNBcBB-1aj-cjxj其中,NB为非基变量的指标集。记为zj-cj=cBB-1xj-cj为检验数,于是有f-f*-jNBzj-cjxj变换后的问题为minfx=f*-jNBzj-cjxjs.t.xB B-1NxNb*xj0(51)其中,f*为基本可行解x*所对应的目标函数值。若基本可行解的所有基变量都取正值,则称它为非退化的; 若有取零值的基变量,则称它为退化的。称所有基本可行解非退化的线性规划为非退化的。对于非退化的线性规划公式(51),有下面结论:(1) 如果所有zj-cj0,则x*为最优解,记为x*。(2) 如果zk-ck0,kNB,且相应的B-1ak0,则无有界最优解。(3) 如果zk-ck0,kNB,且a*k=B-1ak至少有一个正分量,则能找到基本可行解x^,使目标数值下降,有cx^zk-ck=max{zj-cj|j=1,2,,n}确定下标k,取xk为进基变量。(4) 如果zk-ck0,则停止,此时基本可行解x=xBxN=b*0是最优解,目标函数最大值为f=cBb*; 否则,执行下一步。(5) 计算a*k=B-1ak,若a*k0,则停止。此时问题无有界解; 否则执行下一步。(6) 求最小比。b*ra*rk=minb*ia*ik|a*ik0确定下标r,取xBr为离基变量。(7) 以ak代替aBr得到新基,并令xk=b*ra*rk,再返回(2)。在MATLAB中没有提供函数直接实现单纯形算法,下面编写simplefun.m函数实现单纯形算法。代码如下。function [x,f]=simplefunc,A,b%c为线性规划问题%A为其系数矩阵%b为其约束条件%x为最优解%f为最优解处的函数值t=findbAt,: =-At,: ;[m,n]=sizeA; B=A: ,1: m;x=zerosn,1; m1=1: m;whiledetB==0| ~isemptyfindinvB*btp=randpermn;m1=tp1: m;B=A: ,m1;endxB=B\b; xm1=xB;f=cm1*xm1;co=cm1*B\A-c;[z1,z2]=maxco;while z10az=B\A: ,z2;if azdisp''问题无解'',break;elset1=findaz0;[tt1,tt2]=minxBt1.azt1;t3=t1tt2; B: ,t3=A: ,z2;xm1=xB-tt1*az;m1t3=z2; xz2=tt1;f=cm1*xB; xB=xm1;co=cm1B*A-c;[z1,z2]=maxco;endendxm1=xB; f=c*x;【例53】利用单纯形算法计算下列线性规划问题:max fz=x1 2x4 3x5s.t.x1 2x2 3x3=8-x1 2x2 4x4 5x5=3x1 4x2-8x4 3x5 5x6=5xj0j=1,2,,6其实现的MATLAB代码如下。 clear all;c=[1 0 0 2 3 0];b=[8 3 5]'';A=[1 2 3 5 0 0; -1 2 0 4 5 0; 1 4 0 8 3 5];x=simplefunc,A,b运行程序,输出如下。x =01.14291.904800.14290如果在基本可行解中存在有基变量为零,则称为退化的基本可行解。在基本可行解退化时,有可能发生用单纯形法进行无限多次迭代也得不到最优解的死循环。如果目标函数的系数全部大于0,则要注意约束的问题,可用以下自定义编写的函数实现,代码如下。function [x,f]=CmpSimpleMthA,c,b,baseV%约束矩阵A%目标函数系数向量c%约束右端向量b%初始基向量baseV%目标函数取最小值时的自变量值x%目标函数的最小值fsz=sizeAnV=sz2;n=sz1;xx=1: nV;nob=zeros1,1;m=1;if c=0vr=findc~=0,1,''last'';rgv=invA: ,nV-n 1: nV*b;if rgv=0x=zeros1,vr;f=0;elsedisp''不存在最优解!'';x=NaN;f=NaN;return;endendfor i=1: nV%获取非基变量下标ifisemptyfindbaseV==xxi,1nobm=i;m=m 1;else;endendbCon=1;M=0;while bConnB=A: ,nob; %非基变量矩阵ncb=cnob; %非基变量系数B=A: ,baseV; %基变量矩阵cb=cbaseV; %基变量系数xb=invB*b;ff=cb*xb;w=cb*invB;for i=1: lengthnob %判别sigmai=w*nB: ,i-ncbi;end[maxs,ind]=maxsigma; %ind为进基变量下标if maxsminf=cb*xb;vr=findc~=0,1,''last'';for l=1: vrele=findbaseV==l,1;ifisemptyelexl=0;elsexl=xbele;endendbCon=0;elsey=invB*A: ,nobind;if ydisp''不存在最优解!'';x=NaN;f=NaN;return;else%寻找出基变量minb=inf;chagB=0;for j=1: lengthyif yj0bz=xbjyj;if bzminb=bz;chagB=j;endendend%chagB为出基变量下标tmp=baseVchagB; %更新基矩阵和非基矩阵baseVchagB=nobind;nobind=tmp;endendM=M 1;ifM==1000000 %迭代步数限制disp''找不到最优解!'';x=NaN;f=NaN;return;endend【例54】用单纯形法求解下面的线性规划。min f=2x1 x2s.t.-x1-2x242x1 3x212x1,x20将上面的一般形式化为标准形式,有min f=2x1 x2s.t.-x1-2x2-x3=42x1 3x2 x4=12x1,x2,x3,x4,x50其实现的MATLAB代码为 clear all;A=[-1 -2 -1 0; 2 3 0 1];c=[2 1 0 0];b=[4; 12];[x,f]=CmpSimpleMthA,c,b,[2,3]运行程序,输出如下。不存在最优解!x =NaNf =NaN这样得到的结果才是正确的。5.4修正单纯形法在单纯形算法中,每步都需要计算基矩阵B的逆,而逆矩阵的计算是很花时间的。修正单纯形法通过对旧的矩阵B的逆作行变换,来得到新的矩阵B的逆,因此只需要在迭代的初始计算矩阵B的逆,然后通过行变换来求逆,这样就缩短了求解时间。修正单纯形法的思想是计算出下式bsysk=minbiyik|yik0然后对原来的B-1作如下的行变换,设B-1=b~ij,令b~ij=b~ij-yikyskb~sjisb~sj=b~sjysk这样就可得到新的逆矩阵子。在MATLAB中,采用自定义的修正单纯形函数ModiMthd.m实现求解,函数的源代码如下。function [x, minf]=ModiMthdA,c,b,baseV%约束矩阵A%目标函数系数向量c%约束右端向量b%初始基向量baseV%目标函数取最小值时的自变量值x%目标函数的最小值 minfsz=sizeAnV=sz2;n=sz1;xx=1: nV;nob=zeros1,1;m=1;if c=0vr=findc~=0,1,''last'';rgv=invA: ,nV-n 1: nV*b;if rgv=0x=zeros1,vr;minf=0;elsedisp''不存在最优解!'';x=NaN;minf=NaN;return;endendfor i=1: nV%获取非基变量下标ifisemptyfindbaseV==xxi,1nobm=i;m=m 1;else;endendbCon=1;M=0;B=A: ,baseV;invB=invB;while bConnB=A: ,nob; %非基变量矩阵ncb=cnob; %非基变量系数B=A: ,baseV; %基变量矩阵cb=cbaseV; %基变量系数xb=invB*b;f=cb*xb;w=cb*invB;for i=1: lengthnob %判别sigmai=w*nB: ,i-ncbi;end[maxs,ind]=maxsigma; %ind为进基变量下标if maxsminf=cb*xb;vr=findc~=0,1,''last'';for l=1: vrele=findbaseV==l,1;ifisemptyelexl=0;elsexl=xbele;endendbCon=0;elsey=invB*A: ,nobind;if ydisp''不存在最优解!'';x=NaN;minf=NaN;return;else%寻找出基变量minb=inf;chagB=0;for j=1: lengthyif yj0bz=xbjyj;if bzminb=bz;chagB=j;endendend%chagB为出基变量下标tmp=baseVchagB; %更新基矩阵和非基矩阵baseVchagB=nobind;nobind=tmp;for j=1: chagB-1 %基变量矩阵的逆矩阵变换if yj~=0;invBj,: =invBj,: -invBchagB,: *yjychagB;endendfor j=chagB 1: lengthyif yj~=0invBj,: =invBj,: -invBchagB,: *yjychagB;endendinvBchagB,: =invBchagB,: ychagB;endendM=M 1;ifM==1000000 %迭代步数限制disp''找不到最优解!'';x=NaN;minf=NaN;return;endend【例55】用修正单纯形法求解下面的线性规划问题。min f=x1-3x2 x3s.t.2x1-x2 x3=82x1 x22x1 2x210x1,x2,x30将上面一般形式化为标准形式min f=x1-3x2 x3s.t.2x1-x2 x3=82x1 x2-x4=2x1 2x2 x5=10x1,x2,x3,x4,x50得出:A=2-1100210-1010001,b=8210,c=[1-3100]其实现的MATLAB代码为 clear all;A=[2 -1 1 0 0; 2 1 0 -1 0; 1 2 0 0 1];c=[1 -3 1 0 0];b=[8 2 10]'';baseV=[3 4 5];[x,f]=ModiMthdA,c,b,baseV运行程序,输出如下。sz =3 5x =0 513f =-25.5大M法大M法求解线性规划的过程与单纯形法一样,不同的是对线性规划一般形式的处理方式,大M法将线性规划转化成如下问题:minf=cxs.t.Ax=bx0minf=cx MeTys.t.Ax=bx0在MATLAB中,采用自定义的大M法函数simple_Mfun.m实现求解,函数的源代码如下。function [x,f]=simple_Mfunc,A,b,M,N,eps%M为充分大的数%N为引进入人工变量的个数,N应不超过通常等于约束等式的个数%eps为精度%x为最优解%f为最优解处的函数值[m,n]=sizeA;if nargineps=0;endif narginN=0;endif NMerror''N不能超过约束条件的个数m!!!'';elseA=[A,[zerosN,m-N; eyem-N]];c=[c: '',zeros1,m-N];A=[A,eyem];c=[c,M*ones1,m];m1=n m-N 1: n 2*m-N;B=A: ,m1;x=zerosn 2*m-N,1;x=x: ;t=findbAt,: =-At,: ; xB=B\b;xm1=xB; f=c*x;co=cm1*B\A-c; [z1,z2]=maxco;whilez1epsaz=B\A: ,z2;if azx=nan*oneslengthc;break;elset1=findazeps;p=[xB,B\eyesizeB];pp=[];for k=1: lengtht1ppk,: =pt1k,: .azt1k;endtt1=minxBt1.azt1;[tt0,tt2]=min_Mpp;t3=t1tt2; B: ,t3=A: ,z2;xm1=xB-tt1*az;m1t3=z2; xz2=tt1;f=cm1*xB; xB=xm1;co=cm1*B\A-c; [z1,z2]=maxco;endendendif sumxn m-N 1: n 2*m-Nxm1=xB; f=c1: n*x1: n;x=x1: n;elsex=nan*oneslengthc;x=x: ; x=x1: n;end在以上编写的simple_Mfun.m函数过程中,调用到采用字典序最小法编写的min_M.m函数。其代码如下。function [y,k]=min_Mx%x线性规划问题%y的返回值为矩阵x字典序最小行向量%k为y在x中的行数[m,n]=sizex;k=1; y=x1,: ;ifm==1,k=1; y=x1,: ;else[t1,t2]=minx;t3=zerosm,n;for i=1: nt3: ,i=t1i*onesm,1;endt4=sumt3~=x;t5=findt4~=0;k=t2t51; y=xk,: ;end【例56】采用大M法求解如下线性规划问题,看是否出现死循环。min fz=-0.75x4 20x5-0.5x6 6x7s.t.x1 0.25x4-8x5-x6 9x7=0x2 0.5x4-12x5-0.6x6 9x7=0x3 x6=1xj0j=1,2,,7其实现的MATLAB代码如下。 clear all;c=[0 0 0 -0.75 20 -0.5 6];a=[0.25 -8 -1 9; 0.5 -12 -0.5 3; 0 0 1 0];A=[eye3,a];b=[0 0 1]'';lb=zeros7,1;x=linprogc,[],[],A,b,lb;xx=simple_Mfunc,A,b,100000,3;[x,xx]运行程序,输出如下。Optimization terminated.ans =0.75000.75000.0000 00.0000 01.00001.00000.0000 01.00001.00000.0000 -0.0000可见采用simple_Mfun函数与linprog函数计算结果是一样的,证明没有出现死循环。5.6有界单纯形法前面介绍的几种算法都是基于自变量非负的假设,而实际的线性规划问题,其自变量往往只能在某个区间取值,即lbxub。即变量有界单纯形法如下形式的线性规划:min f=cxs.t.Ax=blbxub变量有界单纯形法求解线性规划的算法与单纯形法求解线性规划的算法相同,在MATLAB中都是通过自定义的simplefun.m函数实现。【例57】用有界单纯形法求解下面的线性规划问题:min f=-2x1 x2s.t.-x1 2x242x1 3x212x12,x210把以上的线性规划问题转化为标准形式,得min f=-2x1 x2s.t.-x1 2x2-x3=42x1 3x2-x4=12x1-x5=2x1 x6=10x2-x7=2x2 x8=10x1,x2,x3,x4,x5,x6,x7,x80得出:A=-12-100000230-100001000-100010000100010000-1001000001c=[-21000000]b=[412210210]T其实现的MATLAB代码为 clear all;A=[-1 2 -1 0 0 0 0 0; 2 3 0 -1 0 0 0 0; 1 0 0 0 -1 0 0 0; 1 0 0 0 0 1 0 0; 0 1 0 0 0 0 -1 0; 0 1 0 0 0 0 0 1];b=[4 12 2 10 2 10]'';c=[-2 1 0 0 0 0 0 0];[x,f]=simplefunc,A,b运行程序,输出结果如下。x =1070298053f =-135.7MATLAB函数实现线性规划与一般的线性规划理论一样,在MATLAB中有线性规划的标准型,和前面介绍的一般标准型有类似之处,也有不同之处。在调用MATLAB线性规划函数linprog时,要遵循MATLAB中对标准型的要求。5.7.1MATLAB线性规划标准型线性规划问题的MATLAB标准型为min f=cTxs.t.AxbAeqx=beqlbxub在上述模型中,有一个需要极小化的目标函数f,以及需要满足的约束条件。假设x为n维设计变量,且线性规划问题具有不等式约束m1个、等式约束m2个,那么c、x、lb和ub均为n维列向量,b为m1维列向量,beq为m2维列向量,A为m1n维矩阵,Aeq为m2n维矩阵。需要注意的是如下两点:(1) 在该MATLAB标准型中,目的是对目标函数求极小;(2) MATLAB标准型中的不等式约束形式为。【例58】对于如下线性规划问题:max f=4x1-2x2 x3s.t.2x1-x2 x312-8x1 2x2-2x38-2x1 x3=3x1 x2=7x1,x2,x30要转化为MATLAB标准形式,则需要经过如下几个步骤。(1) 原问题是对目标函数求极大,因此添加负号使问题目标为minf=-4x1 2x2-x3。(2) 原问题中存在的约束条件,因此添加负号使其变为8x1-2x2 2x38。于是不等式组合约束写成矩阵的形式2-118-22x1x2x312-8将等式写成矩阵的形式-201110x1x2x3375.7.2MATLAB实现线性规划在MATLAB中,提供实现线性规划问题的函数为linprog函数,该函数的调用格式如下。x = linprogf,A,b: 求minf''*x在约束条件A.xb下线性规划的最优解。x = linprogf,A,b,Aeq,beq: 等式约束Aeq.x=beq,如果没有不等式约束A.xb,则置A=[],b=[]。x = linprogf,A,b,Aeq,beq,lb,ub: 指定x的范围lbxub,如果没有等式约束Aeq.x=beq,则置Aeq=[],beq=[]。x = linprogf,A,b,Aeq,beq,lb,ub,x0: x0为给定初始值。x = linprogf,A,b,Aeq,beq,lb,ub,x0,options: options为指定的优化参数,如表53所示。表53options优化参数优化参数说明LargeScale如果设置为on,则使用大规模算法; 如果设置为off,则使用中小规模算法Diagnostics打印要极小化的函数的诊断信息Display设置为off,则不显示输出; 为iter,则显示每一次的迭代输出; 为final,则只显示最终结果MaxIter函数所允许的最大迭代次数TolFun函数值的容忍度[x,fval] = linprog: fval为返回目标函数的最优值,即fval=c''x。[x,fval,exitflag] = linprog: exitflag为终止迭代的错误条件,其参数如表54所示。表54exitflag的值及说明exitflag的值说明1表示函数收敛到解x0表示达到了函数最大评价次数或迭代的最大次数-2表示没有找到可行解-3表示所求解的线性规划问题是无界的-4表示在执行算法时遇到了NaN-5表示原问题和对偶问题都是不可行的-7表示搜索方向使得目标函数数值下降的很少[x,fval,exitflag,output] = linprog: output为关于优化的一些信息,其结构及说明如表55所示。表55output结构及说明output结构说明iterations表示算法的迭代次数algorithm表示求解线性规划问题时所用的算法cgiterations表示共轭梯度迭代的次数message表示算法的退出信息firstorderopt一阶最优测量constrviolation最大约束函数[x,fval,exitflag,output,lambda] = linprog: lambda为输出各种约束对应的Lagrange乘子(即为相应的对偶变量值),它是一个结构体变量,其结构及说明如表56所示。表56lambda结构及说明lambda结构说明ineqlin表示不等式约束对应的Lagrange乘子向量eqlin表示等式约束对应的Lagrange乘子向量upper表示上界约束xub对应的Lagrange乘子向量lower表示下界约束xlb对应的Lagrange乘子向量下面通过实例演示linprog函数的用法。【例59】对下列非线性方程进行线性规划问题求解,min fx=-5x1-4x2-6x3s.t.x1-x2 x3203x1 2x2 4x3423x1 2x230x1,x2,x30其实现的MATLAB代码为 clear all;f = [-5; -4; -6]; %输入目标函数系数矩阵A =[1 -11; 324; 320]; %输入不等式约束系数矩阵b = [20; 42; 30]; %输入右端点lb = zeros3,1;[x,fval,exitflag,output,lambda]=linprogf,A,b,[],[],lb,[],[],optimset''Display'',''iter''运行程序,输出如下。Residuals: Primal Dual DualityTotalInfeasInfeasGap RelA*x-bA''*y z-fx''*zError--------------------------------------------Iter0: 1.13e 003 1.87e 001 2.62e 003 1.50e 003Iter1: 1.64e 002 2.05e-015 4.34e 002 2.96e 000Iter2: 2.95e-014 2.00e-015 7.92e 001 5.41e-001Iter3: 5.02e-015 8.24e-015 7.70e 000 9.50e-002Iter4: 7.79e-012 2.19e-015 7.74e-001 9.84e-003Iter5: 4.16e-014 1.56e-015 2.01e-004 2.57e-006Iter6: 1.00e-014 1.78e-015 2.01e-009 2.57e-011Optimization terminated.x = %最优解0.000015.00003.0000fval = %最优值-78.0000exitflag =%算法收敛1output =iterations: 6 %迭代次数algorithm: ''large-scale: interior point'' %应用的是内点算法cgiterations: 0 %没有共轭梯度迭代message: ''Optimization terminated.'' %达到最优解,终止迭代constrviolation: 0 %没有最大约束函数firstorderopt: 5.8704e-010 %一阶最优测量lambda =%Lagrange乘子向量ineqlin: [3x1 double]eqlin: [0x1 double]upper: [3x1 double]lower: [3x1 double]【例510】求解下面4元线性规划问题,max fx=34x1-150x2 150x3-6x4s.t.14x1-60x2-150x3 9x40-12x1 90x2 150x3-3x40x1,x2,x4-5,-5x31原问题中求解的是最大值问题,将其转换为最小化问题:max fx=-34x1 150x2-150x3 6x4s.t.14x1-60x2-150x3 9x4012x1-90x2-150x3 3x40x1,x2,x4-5,-5x31其实现的MATLAB代码为clear all;f = [-34 150 -150 6];A=[14 -60 -150 9; 12 -90 -150 3];b=[0 0]'';Aeq=[]; beq=[];lb=[-5 -5 -5 -5]''; %下界ub=[inf inf 1 inf]''; %上界[x,fval,exitflag,output,lambda]=linprogf,A,b,Aeq,beq,lb,ub,[0 0 0 0]'',optimset''Display'',''iter''运行程序,输出如下。Residuals: Primal DualUpperDuality TotalInfeasInfeasBounds GapRelA*x-b A''*y z-w-f {x} s-ubx''*z s''*w Error------------------------------------------------------Iter0: 9.39e 003 1.43e 002 1.94e 002 6.03e 004 2.77e 001Iter1: 5.90e-012 1.21e 001 0.00e 000 3.50e 003 1.78e 000Iter2: 3.22e-011 5.01e-001 0.00e 000 3.61e 002 3.79e-001Iter3: 1.16e-010 2.24e-001 0.00e 000 3.77e 002 3.62e-001Iter4: 1.35e-011 7.22e-016 0.00e 000 3.37e 001 4.59e-002Iter5: 1.80e-013 1.53e-015 0.00e 000 6.86e-001 9.51e-004Iter6: 7.19e-013 1.78e-015 0.00e 000 2.38e-002 3.30e-005Iter7: 5.86e-013 2.84e-014 8.88e-016 1.19e-006 1.66e-009Optimization terminated.x =-5.0000-0.19471.0000-5.0000fval =-55.4700exitflag = 1output =iterations: 7algorithm: ''large-scale: interior point''cgiterations: 0message: ''Optimization terminated.''constrviolation: 0lambda =ineqlin: [2x1 double]eqlin: [0x1 double]upper: [4x1 double]lower: [4x1 double]5.8线性规划的实际应用本节中的实例将线性规划问题的建模和求解综合在一起,目的是让读者进一步了解典型线性规划问题的建模技巧和MATLAB在线性规划问题求解中的使用方法。5.8.1生产决策问题【例511】某车间生产A、B两种产品,已知生产产品A、B需要用原料分别为3个单位和4个单位,所需的工时分别为5个单位和3个单位,现在可以应用的原料为120个单位,工时为150个单位,每生产一台A和B分别可获得7元和5元,应当如何安排生产A、B的台数,才能使车间获得最大利润?令生产产品A、B分别为x1、x1台,由题可建立如下数学模型:max fx=7x1 5x2s.t.2.5x1 3.5x21205x1 3x2150x1,x20将该数学模型按照MATLAB的要求转换为目标函数最小化,即max fx=-7x1-5x2s.t.2.5x1 3.5x21205x1 3x2150x1,x20其实现的MATLAB代码如下。 clear all;f=[-7,-5]''; %目标函数A=[2.5,3.5; 5,3]; %约束方程中变量的系数b=[120,150]''; %约束方程的系数lb=[0,0]''; %变量下界[x,fval,exitflag,output,lambda]=linprogf,A,b,[],[],lb %线性规则问题求解运行程序,输出如下:Optimization terminated.x =16.500022.5000fval =-228.0000exitflag =1output =iterations: 5algorithm: ''interior-point''cgiterations: 0message: ''Optimization terminated.''constrviolation: 0firstorderopt: 2.3119e-07lambda =ineqlin: [2x1 double]eqlin: [0x1 double]upper: [2x1 double]lower: [2x1 double]由结果可知,生产产品A、B的数量分别为17台、22台或生产产品A、B的数量分别为16台、23台时利润最大,最大利润为228元。5.8.2生产计划安排问题【例512】某工厂计划生产甲、乙两种产品,主要材料有钢材3500kg,铁材1800kg、专用设备能力2800台时,材料与设备能力的消耗定额及单位产品所获利润如表57所示,问如何安排生产,才能使该厂所获利润最大?表57材料与设备能力的消耗及单位产品所获利润产品单位产品消耗定额材料与设备甲件乙件现在材料与设备能力钢材kg853500铁材kg641800设备能力台时452800单位产品的利润元80125首先建立模型,设甲、乙两种产品计划生产量分别为x1、x2(件),总的利润为fx(元)。求变量x1、x2的值为多少时,才能使总利润fx=80x1 125x2最大?依题意可建立数学模型为maxfx=80x1 125x2s.t.8x1 5x235006x1 4x218004x1 5x22800x1,x20x10,x20,x30因为linprog是求极小值问题,所以以上模型可变为maxfx=-80x1-125x2s.t.8x1 5x235006x1 4x218004x1 5x22800x1,x20根据上述模型,其实现的MATLAB代码如下。 clear all;F=[-80,-125];A=[8 5; 6 4; 4 5];b=[3500,1800,2800];lb=[0; 0]; ub=[inf; inf];[x,fval, exitflag, output] = linprogF,A,b,[],[],lb%线性规划问题求解运行程序,输出如下。Optimization terminated.x =0.0000450.0000fval =-5.6250e 004exitflag =1output =iterations: 5algorithm: ''large-scale: interior point''cgiterations: 0message: ''Optimization terminated.''constrviolation: 0firstorderopt: 2.2804e-10当决策变量x=x1,x2=0,450时,规划问题有最优解,此时目标函数的最小值是fval=56250,即当不生产甲产品,只生产乙产品450件时,该厂可获最大利润为56250元。5.8.3配料问题【例513】某种作物在全部生产过程中至少需要32kg氮,磷以24kg为宜,钾不得超过42kg。现有甲、乙、丙、丁4种肥料,各种肥料的单位价格及含氮、磷、钾数量如表58所示。表58各种肥料的单位价格及含氮、磷、钾数量单位: kg各种元素及价格甲乙丙丁氮0.030.300.15磷0.0500.200.10钾0.14000.07价格0.040.150.100.125请问,应如何配合使用这些肥料,使得既能满足作物对氮、磷、钾的需要,又能使施肥成本最低?假设以决策变量x1、x2、x3、x4分别表示甲、乙、丙、丁4种肥料的用量,从而根据表58得到如下的线性方程组minfz=0.04x1 0.15x2 0.1x3 0.125x4s.t.0.03x1 0.3x2 0.15x4320.05x1 0.2x3 0.1x4=240.14x1 0.07x442x1,x2,x3,x40将约束条件化为标准的线性规划形式:minfz=0.04x1 0.15x2 0.1x3 0.125x4s.t.-0.03x1-0.3x2-0.15x4-320.05x1 0.2x3 0.1x4=240.14x1 0.07x442x1,x2,x3,x40其实现的MATLAB代码为 clear all;f=[0.04 0.15 0.1 0.125]''; %目标函数A=[-0.03 -0.3 0 -0.15; 0.14 0 0 0.07];b=[-32 42]'';Aeq=[0.05 0 0.2 0.1];Beq=[24];lb=zeros4,1; %优化变量的下限,无上限约束[x,fval,exitflag]=linprogf,A,b,Aeq,Beq,lb%线性规划问题求解运行程序,输出如下。Optimization terminated.x =76.096561.743663.662574.6268fval =28.0000exitflag =1由以上结果可看出,当甲、乙、丙、丁4种肥料的用量分别为76.0965、61.7436、63.6625、74.6268kg时,即既能满足作物对氮、磷、钾的需要,又能使施肥成本最低,为28元。5.8.4投资选择问题【例514】某投资者有50万元资金可用于长期投资,可供选择的投资品种包括购买国债、公司债券、股票、银行储蓄与投资房地产。各种投资方式的投资期限、年收益率、风险系数、增长潜力的具体参数如表59所示。若投资者希望投资组合的平均年限不超过5年,平均的期望收益率不低于12.5%,风险系数不超过3.5,收益的增长潜力不低于10%。问在满足上述要求的条件下,投资者该如何进行组合投资,使平均年收益率达到最高?表59各种投资方式的投资期限、年收益率、风险系数、增长潜力序号投资方式投资年限年年收益率%风险系数增长潜力%1国债311102公司债券8143163房地产5219304股票3248245短期储蓄160.526长期储蓄4151.54首先,建立目标函数。设决策变量为x1、x2、x3、x4、x5、x6,其中xi为第i种投资方式在总投资中占的比例,由于决策的目标是使投资组合的平均年收益率最高,因此目标函数为max fx=11x1 14x2 21x3 24x4 6x5 15x6再根据题意建立约束条件s.t.3x1 8x2 5x3 3x4 x5 4x6511x1 14x2 21x3 24x4 6x5 15x612.5x1 3x2 9x3 8x4 0.5x5 1.5x63.516x2 30x3 24x4 2x5 4x610x1 x2 x3 x4 x5 x6=1x1,x2,x3,x4,x5,x60即其数学模型为max fx=11x1 14x2 21x3 24x4 6x5 15x6s.t.3x1 8x2 5x3 3x4 x5 4x6511x1 14x2 21x3 24x4 6x5 15x612.5x1 3x2 9x3 8x4 0.5x5 1.5x63.516x2 30x3 24x4 2x5 4x610x1 x2 x3 x4 x5 x6=1x1,x2,x3,x4,x5,x60根据linprog函数的要求将数学模型改为min fx=-11x1-14x2-21x3-24x4-6x5-15x6s.t.3x1 8x2 5x3 3x4 x5 4x65-11x1-14x2-21x3-24x4-6x5-15x6-12.5x1 3x2 9x3 8x4 0.5x5 1.5x63.5-16x2-30x3-24x4-2x5-4x6-10x1 x2 x3 x4 x5 x6=1x1,x2,x3,x4,x5,x60其实现的MATLAB代码如下。 clear all;c=[-11 -14 -21 -24 -6 -15];A=[3 8 5 3 1 4; -11 -14 -21 -24 -6 -15; 1 3 9 8 0.5 1.5; 0 -16 -30 -24 -2 -4];b=[5 -12.5 3.5 -10]'';lb=zeros6,1;Aeq=ones1,6;beq=1;[x,fval]=linprogc,A,b,Aeq,beq,lb运行程序,输出如下。Optimization terminated.x =0.00000.00000.00000.30770.00000.6923fval =-17.7692运行结果表明,投资组合选择的决策是长期储蓄占投资总额的69.23%,股票投资占总额的30.77%。其年收益为17.7692万元。5.8.5平衡指派问题在现实生活中,有各种性质的指派问题。例如,有若干项工作需要分配给若干人(或部分)来完成,有若干项合同需要选择若干个投标者来承包,有若干班级需要安排在各教案上课等。诸如此类的问题,它们的基本要求是在满足特定的指派要求条件下,使指派方案的总体效果最佳。由于指派问题的多样性,有必要定义指派问题的标准形式。指派问题的标准形式(以人和事为例)是: 有n个人和n件事,已知第i个人做第j件事的费用为ciji,j=1,2,,n,要求确定人和事之间的一一对应的指派方案,使完成这n件事的总费用最少。为了建立标准指派问题的数学模型,引入n2个01变量:xij=0如果不能派第i人做第j件事件,i,j=1,2,,n1如果指派第i人做第j件事件,i,j=1,2,,n这样,问题的数学模型可写为:min z=ni=1nj=1cijxijs.t.ni=1xij=1j=1,2,,nnj=1xij=1i=1,2,,nxij=0,1j=1,2,,n其中,第一个等式约束条件表示每件事必由且只有一个人去做,第二个等式约束条件表示每个人必做且只做一件事。指派问题是产量(ai)、销量(bj)相等,且ai=bj=1i,j=1,2,,n的运输问题。有时也称cij为第i个人完成第j件工作所需的资源数,称之为效率系数(或价值系数),并称矩阵C=cijnn=c11c12c1nc21c22c2ncn1cn2cnn为效率矩阵(或价值系数矩阵),称决策变量xij排成的nn矩阵X=xijnn=x11x12x1nx21x22x2nxn1xn2xnn为决策变量矩阵。其特征是有n个1,其他都是0。这n个1位于不同行、不同列。每一种情况为指派问题的一个可行解,共n!个解。除了上述一个人只能做一件事之外,指派问题还有其他的类型。如若某个人可做几件事,则可将该人化作相同的几个人来接近指派。这几个人做同一件事的费用系数当然一样。对于人数和事数不等的指派问题,如果人数小于事数,则添一些虚拟的人,此时这些虚拟的人做各件事的费用系数取为0,理解为这些费用实际上不会发生。如果人数大于事数,则添一些虚拟的事,此时这些虚拟的事被各个人做的费用系数同样也取为0。如果某事不能由某人去做,可将此人做此事的费用取作足够大的M。计算指派问题,首先考虑扩展的指派问题模型:min z=ni=1nj=1cijxijs.t.ni=1xij=1j=1,2,,nnj=1xij=1i=1,2,,n,ai=bj=10xij1i,j=1,2,,n【例515】4个工人被分派做4项工作,规定每人只能做一项工作,每项工作只能一个人做,现设每个人做每项工作所消耗的时间如表510所示,求总耗时最少的分派方案。表510工作时间表单位: 小时工人名工作1工作2工作3工作4工人115182124工人219232218工人326171619工人419212317其实现的MATLAB代码为 clear all;A=[15 18 21 24; 19 23 22 18; 26 17 16 19; 19 21 23 17];a=A''; %效率函数f=a: ; %目标函数o=ones1,4;z=zeros1,4;y=eye4;aeq=[o z z z; z o z z; z z o z; z z z o];aeq=[aeq; y y y y];beq=ones8,1;lb=zeros16,1;[x,fval,exitflag,output,lambda]=linprogf,[],[],aeq,beq,lbxx=reshapex,4,4;xx=xx''; %创建指派方阵Optimization terminated.x =0.48340.51660.00000.00000.51660.00000.00000.48340.00000.00001.00000.00000.00000.48340.00000.5166fval = 70.0000exitflag = 1output =iterations: 6algorithm: ''large-scale: interior point''cgiterations: 0message: ''Optimization terminated.''constrviolation: 1.1102e-015firstorderopt: 1.6566e-008lambda =ineqlin: [0x1 double]eqlin: [8x1 double]upper: [16x1 double]lower: [16x1 double] xx0=roundxx%对指派方阵取整,并求出最优值xx0 =0 1 0 01 0 0 00 0 1 00 0 0 1最终结果为: 把工作1给工人2,工作2给工人1,工作3给工人3,工作4给工人4,且总的消耗时间为70小时。5.8.6人员安排问题【例516】某昼夜服务的公共交通系统每天各时段(每4小时为一个时段)所需的值班人数如表511所示。这些值班人员在某一时段开始上班后要连续工作8小时(包括轮流用膳时间)。问该公交系统至少需要多少名工作人员才能满足值班的需要?表511值班安排表班次时间段所需人数个106: 0010: 0056210: 0014: 0072314: 0018: 0056418: 0022: 0048522: 0002: 0024602: 0006: 0030首先,设xi为第i个时段开始上班人员数i=1,2,,6。可得上述问题的数学模型为min z=x1 x2 x3 x4 x5 x6s.t.x6 x156x1 x272x2 x356x3 x448x4 x524x5 x630即利用linprog函数求解以上数学模型,代码为 clear all;f=ones6,1;z2=zeros1,2;z3=zeros1,3;z4=zeros1,4;a=[1,z4,1; 1,1,z4; 0,1,1,z3];a=[a; z2,1,1,z2; z3,1,1,0; z4,1,1];a=-a;lb=zeros6,1;b=-[56 72 56 48 24 30]'';[x,fval exitflag,output]=linprogf,a,b,[],[],lb,[]运行程序,输出如下。Optimization terminated.x =40.924231.075829.613118.386911.485918.5141fval =150.0000exitflag =1output =iterations: 6algorithm: ''large-scale: interior point''cgiterations: 0message: ''Optimization terminated.''constrviolation: 5.1159e-13firstorderopt: 1.0374e-12决策变量x为人数,应考虑整数解,为此作 p=roundx''p =413130181119在命令行窗口中输入如下运算: a*p''ans =-[60 72 61 48 29 30]''可以得到,a*p''b,又p''lb,这表明p为可行解。又有p*f= fval =150,由此可见,p=[413130181119]也是最优解。这样得到该公司在6个时段的最优解安排为:1时段安排41人。2时段安排31人。3时段安排30人。4时段安排18人。5时段安排11人。6时段安排19人。总计安排人员150人。5.8.7运输问题回顾例51,这是一个相当典型的线性规划问题。关于运输问题,在线性规划领域有很多讨论,在此先给出例51的MATLAB求解,然后给出运输问题的一般形式和一些结论。【例517】最后给出的线性规划问题的数学模型为min f=10x11 12x12 9x13 8x21 11x22 13x23s.t.x11 x12 x13=35x21 x22 x23=55x11 x21=26x12 x22=38x13 x23=26xij0i=1,2; j=1,2,3该问题符合MATLAB的标准型,其实现的代码为 clear all;%目标函数c=[10 12 9 8 11 13]'';%线性等式约束Aeq=[1 1 1 0 0 0; 0 0 0 1 1 1; 1 0 0 1 0 0; 0 1 0 0 1 0; 0 0 1 0 0 1];beq=[35 55 26 38 26];%设置变量的边界约束,无上界约束lb=[0 0 0 0 0 0]'';ub=[inf inf inf inf inf inf]'';%求最优解x和目标函数值fval,由于无等式约束,因此设置A=[],b=[][x,fval,exitflag]=linprogc,[],[],Aeq,beq,lb,ub运行程序,输出如下。Optimization terminated.x = %最优解向量x0.00009.000026.000026.000029.00000.0000fval = %最优解向量x处的目标函数值869.0000exitflag =1将上述运输问题推广到一般形式,可描述为设某种物资有m个产地A1,A2,,Am和n个销地B1,B2,,Bn,各产地的产量分别为a1,a2,,am,各销地的需求量分别为b1,b2,bn。且由Ai到Bj的单位物资的运价为cij,问应该怎样调运才能使各产地调出的物资总量满足各销地的需求并使得总运费最省?为了建立该线性规划问题的数学模型,引入变量xij,其取值为由产地Ai运往销地Bj的该商品数量,则产、销地之间的运量和单位物资的运价可以由表512所示。表512产地和销地之间的运量和单位运费销地运量和运价产地B1B2Bn产量A1x11c11x12c12x1nc1na1A2x21c21x22c22x2nc2na2Amxm1cm1xm2cm2xmncmnam销量b1b2bn根据与例51一样的分析过程,可得该优化问题的数学描述为minmi=1mj=1cijxijs.t.nj=1xij=aii=1,2,,mmj=1xij=bjj=1,2,,nxij0i=1,2,,m; j=1,2,,n显然是一个线性规划问题,进一步地,对产销平衡的运输问题,由于有以下关系式存在:mj=1bj=mi=1mj=1xij=mj=1mi=1xij=mi=1ai可以证明上面两式有可行解的充分必要条件,该结论的一个推论是任何产销地平衡问题都有最优解。运输问题的约束条件的系数矩阵相当特殊,可用比较简单的计算方法,习惯上称为表上作业法。该方法是单纯形法在求解运输问题时的一种简化方法,是由康托洛维奇和希奇柯克两人独立提出的,在此不展开介绍。5.8.8投资收益与风险问题【例518】市场上有n种资产(股票、债券)Sii=1,2,,n供投资者选择,某公司有的相当大的数额为M一笔资金可用作这一时期的投资。公司财务分析人员对Si种资产进行评估,估算在这一时期内购买Si有平均收益率为ri,并预测出购买Si损失率为qi,考虑到投资越分散,总的风险就越小,公司已确定用这笔资金购买若干种资产时,总体风险可用所投资的Si中的最大一个风险来度量。购买Si要付交易费,费率为pi,并且当购买额不超过给定值ui时,交易费按购买ui计算(不买当然无须付费)。另外,假定同期银行存款利率是r0,且既无交易费又无风险(r0=5%)。已知n=4时的相关数据如表513所示。表513不同资产的收益率、损失率、交易费率与限额Siri%qi%pi%uiS1262.411.2S2221.62199S3245.44.653S3252.76.438试给该公司设计一种投资组合方案,即用给定的资金M,有选择地购买若干种资产或存银行生息,使净收益尽可能大,而总体风险尽可能小。可建立一个确定投资比例向量模型,使资产组合的净收益尽可能大,而总体风险尽可能小。建立模型,设x0、x1、x2、x3、x4分别是银行存款与投资于S0、S1、S2、S3、S4的投资比例系数,由于银行存款既无交易费又没有风险故p0=0,q0=0,总体风险可用所投资的Si中最大的一个风险来设定,于是投资组合总体风险为:F=max0i4{xiqi}由于题设给出M为相当大的一笔资金,为了简化模型,可认为该公司购买每一项资产都超过给定的定值ui,于是资产组合的平均收益率为:R=4i=0xiri-pi为了使得平均收益率尽可能大,而总体风险尽可能小,可采取固定总体风险的一个上界q,使得总体收益取得最大,为此建立如下的线性规划模型:max R=4i=0xiri-pis.t.x0 x1 x2 x3 x4=1xiqiqi=0,1,2,3,4x0,x1,x2,x3,x40总体风险的上界为[0 3],取步长为0.01,计算301种不同风险时的总体收益的最大值及相应的投资比例系数,并给出投资方案的净收益率与风险损失率的关系图。其实现的MATLAB代码为 clear all;c=[-25 -20 -19.4 -18.4 -5];A=[2.4 0 0 0 0; 0 1.6 0 0 0; 0 0 5.4 0 0; 0 0 0 2.7 0; 0 0 0 0 0];A1=[1 1 1 1 1];b1=1;lb=[0 0 0 0 0]'';t=0: 0.01: 3;b=tones5,1,: ;for k=1: 301;[x: ,k,Yk]=linprogc,A,b: ,k,A1,b1,lb;endplott,-Y; %绘出投资方案的净收益率与风险损失率的关系图axis[0 4 0 30]%拟合净收益率与风险损失率的关系h1=polyfitt1: 62,-Y1: 62,1h2=polyfitt62: 251,-Y62: 251,1h3=polyfitt251: 301,-Y251: 301,1%输出26种风险时各种资产的投资比例系数与收益矩阵tz=[t1: 10: 251'' [x: ,1: 10: 251'' -Y1: 10: 251'']]运行程序,输出如下,效果如图51所示。h1 =25.33805.0000h2 =2.1718 19.8189h3 =0.0000 25.0000tz =00.00000.00000.00000.00001.00005.00000.10000.04170.06250.01850.03700.84037.53380.20000.08330.12500.03700.07410.6806 10.06760.30000.12500.18750.05560.11110.5208 12.60140.40000.16670.25000.07410.14810.3611 15.13520.50000.20830.31250.09260.18520.2014 17.66900.60000.25000.37500.11110.22220.0417 20.20280.70000.29170.43750.12960.14120.0000 21.15460.80000.33330.50000.14810.01850.0000 21.54810.90000.37500.56250.06250.00000.0000 21.83751.00000.41670.58330.00000.00000.0000 22.08331.10000.45830.54170.00000.00000.0000 22.29171.20000.50000.50000.00000.00000.0000 22.50001.30000.54170.45830.00000.00000.0000 22.70831.40000.58330.41670.00000.00000.0000 22.91671.50000.62500.37500.00000.00000.0000 23.12501.60000.66670.33330.00000.00000.0000 23.33331.70000.70830.29170.00000.00000.0000 23.54171.80000.75000.25000.00000.00000.0000 23.75001.90000.79170.20830.00000.00000.0000 23.95832.00000.83330.16670.00000.00000.0000 24.16672.10000.87500.12500.00000.00000.0000 24.37502.20000.91670.08330.00000.00000.0000 24.58332.30000.95830.04170.00000.00000.0000 24.79172.40001.00000.00000.00000.00000.0000 25.00002.50001.00000.00000.00000.00000.0000 25.0000图51投资方案的净收益率与风险损失率
|