运筹学实验报告

.运筹学实验报告专业：班级：姓名：学号:指导教师：数学与应用数学专业2015-12-18实验目录一、实验目的 (3)二、实验要求 (3)三、实验内容 (3)1、线性规划 (3)2、整数规划 (6)3、非线性规划 (13)4、动态规划 (115)5、排队论 (19)四、需用仪器设备 (26)五、MATLAB优化工具箱使用方法简介 (26)六、LINGO优化软件简介 (26)七、实验总结 (27)一、实验目的1、会利用适当的方法建立相关实际问题的数学模型；2、会用数学规划思想及方法解决实际问题；3、会用排队论思想及方法解决实际问题；4、会用决策论思想及方法解决实际问题；5、掌握MATLAB、LINGO等数学软件的应用；二、实验要求1、七人一组每人至少完成一项实验内容；2、每组上交一份实验报告；3、每人进行1~2分钟实验演示；4、实验成绩比例：出勤：40%课堂提问：20%实验报告：30%实验演示：10%。

三、实验内容1、线性规划例运筹学74页14题Min z=-2x1-x22x1+5x2≤60x1+x2≤183x1+x2≤44X2≤10X1,x2≥0用matlab运行后得到以下结果:the program is with the linear programmingPlease input the constraints number of the linear programming m=6m =6Please input the variant number of the linear programming n=2n =2Please input cost array of the objective function c(n)_T=[-2,-1]'c =-2-1Please input the coefficient matrix of the constraints A(m,n)=[2,5;1,1;3,1;0,1;-1,0;0,-1]A =2 51 13 10 1-1 00 -1Please input the resource array of the program b(m)_T=[60,18,44,10,0,0]'b =60184410Optimization terminated.The optimization solution of the programming is:x =13.00005.0000The optimization value of the programming is:opt_value =-31.0000LINDO程序在命令窗口键入以下内容:max -2x-ysubject to2x+5y<=60x+y<=183x+y<=44y<=10end按solve键在reports window出现:Global optimal solution found.Objective value: 0.000000Total solver iterations: 0Variable Value Reduced CostX 0.000000 2.000000Y 0.000000 1.000000Row Slack or Surplus Dual Price1 0.000000 1.0000002 60.00000 0.0000003 18.00000 0.0000004 44.00000 0.0000005 10.00000 0.000000 2、整数规划课本第二章79页1题Max z=100x1+180x2+70x3s.t. 40x1+50x2+60x3≤100003 x1+6x2+ 2x3≤600x1≤130X2≤80x3≤200x1 x2 x3≥0程序运行及结果：biprogramthe program is with the binary linear programmingPlease input the constraints number of the programming m=5m=5Please input the variant number of the programming n=5n=5Please input cost array of the objective function c(n)_T=[100,180,70]'c=10018070Please input the coefficient matrix of the constraints A(m,n)=[40,50,60;3, 6,2;1,0,0;0,1,0;0,0,1]A=40506036 210001000 1Please input the resource array of the program b(m)_T=[10000;600;130; 80;200]b=1000060013080200Optimization terminated.The optimization solution of the programming is:x=The optimization value of the programming is:opt_value=程序名：intprogram b程序说明：%the programm is with the integer linear programming use branch and bound method!%这个程序是用分支定界法解决整数规划问题%please input the parameters in the main function in the comman d winows%请在命令窗口输入这个主要定义函数的参数function[x,f]=ILp(c,A,b,vlb,vub,x0,neqcstr,pre)%min f=c'*x,s.t.A*x<=b,vlb<=x<=vub%f的最小值等于c的转置乘以x，A乘以x小于等于b,x大于等于vlb小于等于vub %the vectors of x is required as integers as whole%x是整个的整数需要%x0is the initialization,'[]'is also ok%x0是初始值,"[]"也可以是。

%neqcstr is the number of equational constraints,when0can be d elete%neqcstr是平均约束条件的数目，当0能删除时%pre is the concise rate%pre是简明率%x is the integer optimization and f is the optimal value%x是整数规划，f是最优值%%%%%%%%%%%%%%%%if nargin<8,pre=0;%nargin is the factually input variants n umber（这个参数是实际输入的变量个数）if nargin<7,neqcstr=0;if nargin<6,x0=[];if nargin<5,vub=[];if nargin<4,vlb=[];endendendendend%%%%%%%%%%%%%%%%%%%set to column vectors%建立列向量x0=x0(:);c=c(:);b=b(:);vlb=vlb(:);vub=vub(:);mm=1;j=1;nvars=length(c');%number of variants（变量的个数）fvub=inf;xall=[];fall=[];x_f_b=[];[xtemp,ztemp,how]=lp(c,A,b,vlb,vub,x0,neqcstr,-1);ftemp=c'*xtemp;%%%%%%%%%%%%%%%%%%%%%%%if strcmp(how,'ok')%compare between how and ok（how和ok之间的比较）temp0=round(xtemp);%临时变量四舍五入temp1=floor(xtemp);%取比其小的整数temp2=find(abs(xtemp-temp0)>pre);mtemp=length(temp2);if~isempty(temp2)x_f_b=[xtemp;ftemp;vlb;vub];while j<=mmi=1;while i<=mtemp%%%%%%%%%%%%%%%%%%%%%if x_f_b(nvars+1,j)<=fvubvlbl=x_f_b(nvars+2:2*nvars+1,j);vubl=x_f_b(2*nvars+2:3*nvars+1,j);vubl(temp2(i))=temp1(temp2(i));[xtemp,z,how]=lp(c,[A;c'],[b;fvub],vlbl,vubl,x0,neqcstr, -1);ftemp=c'*xtemp;if strcmp(how,'ok')templ0=round(xtemp);templ1=floor(xtemp);templ2=find(abs(xtemp-templ0)>pre);if isempty(templ2)xall=[xall,xtemp];fall=[fall,ftemp];fvub=min([fvub,fall]);elseif ftemp<=fvubx_f_b=[x_f_b,[xtemp;ftemp;vlbl;vubl]];endendend%%%%%%%%%%%%%%%%%%if x_f_b(nvars+1,j)<=fvubvlbr=x_f_b(nvars+2:2*nvars+1,j);vlbr(temp2(i))=temp1(temp2(i))+1;vubr=x_f_b(2*nvars+2:3*nvars+1,j);[xtemp,z,how]=lp(c,[A;c'],[b;fvub],vlbr,vubr,x0,neqc str,-1);ftemp=c'*xtemp;if strcmp(how,'ok')tempr0=round(xtemp);tempr1=floor(xtemp);tempr2=find(abs(xtemp-tempr0)>pre);if isempty(tempr2)xall=[xall,xtemp];fall=[fall,ftemp];fvub=min([fvub,fall]);elseif ftemp<=fvubx_f_b=[x_f_b,[xtemp;ftemp;vlbr;vu br]];endendend%%%%%%%%%%%%%%%%%%%%%i=i+1;end%the second whilexint=x_f_b(1:nvars,:);[m,mm]=size(xint);j=j+1;if j>mmbreakend%the end because the break（因为中断而结束）temp0=round(xint(:,j));temp1=floor(xint(:,j));temp2=find(abs(xint(:,j)-temp0)>pre);mtemp=length(temp2);end%the end of while（结束当前）else%correspond the second if（符合第一个如果）x=xtemp;f=ftemp;end%the end of second if（第二个如果的结束）%%%%%%%%%%%%%%%%%%5if~isempty(fall)fmin=min(fall);nmin=find(fall==fmin);x=xall(:,nmin);f=fmin;endelse%correspond the first if（符合第一个如果）x=nan*ones(1,nvars);endLINDO程序例99页第6题第二问在命令窗口键入以下内容:max -11x1-4x2st-x1+2x2<=45x1+2x2<=162x1-x2<=4endgin x1gin x2按solve键在reports window出现:Global optimal solution found.Objective value: 0.000000Extended solver steps: 0Total solver iterations: 0Variable Value Reduced CostX1 0.000000 11.00000X2 0.000000 4.000000Row Slack or Surplus Dual Price1 0.000000 1.0000002 4.000000 0.0000003 16.00000 0.0000004 4.000000 0.0000003、非线性规划程序名：unpfun1函数unpfun1函数执行实例：（课本第四章152页16题）min4x1+6x2-2x1x2-2x^2,取初始点x0=（1；1）在命令窗口键入以下内容:f = 4*x(1)-2*x(1)*x(2) - x(2)^2-6*x(2);[x,fval]=fminunc(unpfun1,x0)（调用无约束线性规划函数）按运行按钮在solution report 窗口得到以下结果:Warning: Gradient must be provided for trust-region method;using line-search method instead.> In fminunc at 265Optimization terminated: relative infinity-norm of gradient less than options.TolFun.x =1.0e-006 *0.2541 -0.2029fval =1.3173e-0134、动态规划程序名：dynamic; dynfun1_1, dynfun1_2, dynfun1_3;例180页第一题程序说明：dynamic程序：% the programm is with the dynamic programming use the recurisive method for the last to first% this is the main function of the methodfunction[p_opt,fval,u]=dynprog(x,DecisFun,ObjFun,TransFun)% the function is to solve the dynamic example in the textbook% x is the situation variant and its column number represent the stage situation% subfunction DecisFun(k,x) is to solve the decision variant of k stage variant x% subfunction ObjFun(k,x,u) is to stage index function% subfunction TransFun(k,x,u) is the stage transformation function,u is the corresponding decision variant% p_opt has four output,the first is the number of the stage,the second is the optimal road of decision% the third is the optimal stategies of the decision ,the forth is the index function group.% fval is a column vector,is to represent the optimal value correspend to the initial stage is x%k=length(x(1,:));f_opt=nan*ones(size(x));d_opt=f_opt;t_vubm=inf*ones(size(x));x_isnan=~isnan(x);t_vub=inf;%%%%%%%%%%%%%%%%%% to caculate the teminate valuestmp1=find(x_isnan(:,k));tmp2=length(tmp1);for i=1:tmp2u=feval(DecisFun,k,x(i,k));tmp3=length(u);for j=1:tmp3tmp=feval(ObjFun,k,x(tmp1(i),k),u(j));if tmp<=t_vubf_opt(i,k)=tmp;d_opt(i,k)=u(j);t_vub=tmp;endendend%%%%%%%%%%%%%%%%% recurisivefor ii=k-1:-1:1tmp10=find(x_isnan(:,ii));tmp20=length(tmp10);for i=1:tmp20u=feval(DecisFun,ii,x(i,ii));tmp30=length(u);for j=1:tmp30tmp00=feval(ObjFun,ii,x(tmp10(i),ii),u(j));tmp40=feval(TransFun,ii,x(tmp10(i),ii),u(j));tmp50=x(:,ii+1)-tmp40;tmp60=find(tmp50==0);if ~isempty(tmp60),if nargin<5tmp00=tmp00+f_opt(tmp60(1),ii+1);if tmp00<=t_vubm(i,ii)f_opt(i,ii)=tmp00;d_opt(i,ii)=u(j);t_vubm(i,ii)=tmp00;endendendendendendfval=f_opt(tmp1,1); %%%%%%%%%%%%%%%%%%% to write the index and parameter or resultp_opt=[];tmpx=[];tmpd=[];tmpf=[];tmp0=find(x_isnan(:,1));tmp01=length(tmp0);for i=1:tmp01tmpd(i)=d_opt(tmp0(i),1);tmpx(i)=x(tmp0(i),1);tmpf(i)=feval(ObjFun,1,tmpx(i),tmpd(i));p_opt(k*(i-1)+1,[1,2,3,4])=[1,tmpx(i),tmpd(i),tmpf(i)]; for ii=2:ktmpx(i)=feval(TransFun,ii-1,tmpx(i),tmpd(i));tmp1=x(:,ii)-tmpx(i);tmp2=find(tmp1==0);if ~isempty(tmp2)tmpd(i)=d_opt(tmp2(1),ii);endtmpf(i)=feval(ObjFun,ii,tmpx(i),tmpd(i));p_opt(k*(i-1)+ii,[1,2,3,4])=[ii,tmpx(i),tmpd(i),tmpf(i)]; endenddynfun2_1程序：function u=dynfun2_1(k,x)if x==1u=[2,3,4];elseif x==2u=[5;6];elseif x==3u=[5;6];elseif x==4u=[5;6];elseif x==5u=[7;8;9];elseif x==6u=[7;8;9];elseif(x==7)|(x==8)|(x==9)u=10;elseif x==10u=10;enddynfun2_2程序：function v=dynfun2_2(k,x,u)tt=[3;2;1;4;3;1;3;3;5;2;5;1;1;4;2;3;1;5];tmp=[x==1&u==2,x==1&u==3,x==1&u==4,x==2&u==5,x==2&u==6,...x==3&u==5,x==3&u==6,x==4&u==5,x==4&u==6,...x==5&u==7,x==5&u==8,x==5&u==9,x==6&u==7,x==6&u==8,x==6&u==9,x==&u== 10,x==8&u==10,x==9&u==10];v=tmp*tt;dynfun2_3程序：function y=dynfun2_3(k,x,u)y=u;在命令窗口的程序执行过程和结果如下:>>x=nan*ones(3,5);>>x>(1,1)=1;>>x(1:3,2)=(2:4)';>>x(1:2,3)=[5,6]';>>x(1:3,4)=(7:9)';>>x(1,5)=10;>>[p,f]=dynamic(x,'dynfun2_1','dynfun2_2','dynfun2_3') 得到以下结果:p =1 1 3 22 3 5 13 5 7 24 7 10 35 10 10 0f =8该结果表明最短路线为1 3 5 7 10最短路程为85排队论程序：queueM/M/1/k系统课本270页例8.2.5程序及结果如下：queuethe program is with queueing theoryPlease input the system pattern:M/M/1/inf=1,M/M/1/k=2,M/M/c/inf=3,M/M/c/m/m=4! Pattern=2Pt =2（该数字表示选择M/M/1/k系统）Please input the average arrival number in unit time lapta=10lapta =10Please input the average service number in unit time mu=30mu =30Please input the parameter k=2k =2输出结果如下:The service intensity(Untruth) of the system is:ru =0.3333The average queue length is:L =0.3846（平均队长）The average waiting length is:Lq =0.0769（平均等待队长）The lost possibility is:p_k =0.0769The average lost customer number in unit time is:lapta_L =0.7692The truth average arrival customer number in unit time is:lapta_e =9.2308The average delay time is:W =0.0417（平均逗留时间）The average waiting time is:Wq =0.0083（平均等待时间）The average service intensity(efficent and truth) of the system is:ru_e =0.3077若按照甲方案，则如下结果：the program is with queueing theoryPlease input the system pattern:M/M/1/inf=1,M/M/1/k=2,M/M/c/inf=3,M/M/c/m/m=4! Pattern=2Pt =2Please input the average arrival number in unit time lapta=10lapta =10Please input the average service number in unit time mu=30mu =30Please input the parameter k=3k =3输出结果如下:The service intensity(Untruth) of the system is:ru =0.25The average queue length is:L =0.3176The average waiting length is:Lq =0.0706The lost possibility is:p_k =0.0118The average lost customer number in unit time is:lapta_L =0.1176The truth average arrival customer number in unit time is: lapta_e =9.8824The average delay time is:W =0.0321The average waiting time is:Wq =0.0071The average service intensity(efficent and truth) of the system is:ru_e =0.2471the program is with queueing theoryPlease input the system pattern:M/M/1/inf=1,M/M/1/k=2,M/M/c/inf=3,M/M/c/m/m=4! Pattern=2Pt =2Please input the average arrival number in unit time lapta=10lapta =10Please input the average service number in unit time mu=40mu =40Please input the parameter k=2k =2输出结果如下:The service intensity(Untruth) of the system is:ru =0.25The average queue length is:L =0.2857The average waiting length is:Lq =0.0476The lost possibility is:p_k =0.0476The average lost customer number in unit time is:lapta_L =0.4762The truth average arrival customer number in unit time is:lapta_e =9.5238The average delay time is:W =0.03The average waiting time is:Wq =0.005The average service intensity(efficent and truth) of the system is:ru_e =0.2381the program is with queueing theoryPlease input the system pattern:M/M/1/inf=1,M/M/1/k=2,M/M/c/inf=3,M/M/c/m/m=4! Pattern=2Pt =2Please input the average arrival number in unit time lapta=30lapta =30Please input the average service number in unit time mu=30mu =30Please input the parameter k=2k =2输出结果如下:The service intensity(Untruth) of the system is:ru =1The average queue length is:L =1The average waiting length is:Lq =0.3333The lost possibility is:p_k =0.3333The average lost customer number in unit time is:lapta_L =10The truth average arrival customer number in unit time is: lapta_e =20The average delay time is:W =0.05The average waiting time is:Wq =0.0167The average service intensity(efficent and truth) of the system is:ru_e =0.6667四、需用仪器设备PC i5、windows XP、MATLABR2007a、LINGO11五、MATLAB优化工具箱使用方法简介MATLAB优化工具箱具有强大的科学计算能力,在工程设计领域得到了广泛的应用.简要介绍了MATLAB优化工具箱,通过对MATtAB优化工具箱中fmincon函数的语法进行分析,提出了结构优化设计的通用求解方法.首先,合理设置优化目标函数和约束条件.然后,使用MATLAB优化工具箱进行编程计算.结果显示,与其他方法相比,使用MATLAB优化工具箱进行优化,不仅可以提高计算精度,而且可以减少计算时间.因此,在结构工程设计领域有较强的实际应用价值.六、LINGO优化软件简介lingo是用来求解线性和非线性优化问题的简易工具。