递归程序设计与优化模型及其求解.doc

1 递归程序设计递归算法实现的程序简洁、易读懂。可以用来解决一些复杂的问题。在使用Matlab软件进行递归程序设计当中，要特别注意的是：程序中一定要有返回语句(return)，以便使程序在适当的时候终止递归调用，否则会出现堆栈溢出，甚至导致程序中断退出。为防止程序异常终止，可以添加异常处理代码(trycatchend)。递归程序设计的第一步，要分析问题，建立问题求解的递推关系式或算法。然后再进行程序设计实现。1.1 计算阶乘从n个数中取n个数的数的全排列记为，且有。因此建立其求阶乘的递推关系式：，结束条件：。计算阶乘的Matlab函数如下：function r=factorial(n)%计算n的阶乘if n=1 r=1; return endr=n*factorial(n-1);例子：再命令行输入factorial(3)则返回61.2 组合数学中的Pascal公式从n个数中取r个的组合数记为，且有，。组合数学中又这样的Pascal公式：。因此可以编写如下的函数求组合数：function num=getcom(n,r)if r=0 | (n=r) %结束条件 num=1; returnend if r=1,%2004-12-7add num=n; returnendnum=getcom(n-1,r)+getcom(n-1,r-1);1.3 汉诺塔问题1.3.1 “Hanoi塔”问题有3根柱子：A,B,C，现有n个大小不一的圆盘依半径的大小，从下而上套在柱子A上，最大的圆盘放在柱子A的最下面。现要将所有的圆盘从柱子A移动到C柱子上，每次只允许从一根柱子转移到另一根柱子上，且在转移过程中不允许出现大圆盘放在小圆盘上。B盘为可以利用的柱子，每次只允许移动一个盘子，请问要转移多少次才能将柱子A上的圆盘全部转移柱子C上？“Hanoi塔”是组合数学中的著名问题之一。1.3.2 问题求解主程序调用：global nmovenmove=0;hanta(A,B,C,3)nmove说明：上面的程序调用表示有3根柱子：A ,B,C，现有3个盘子在A上，要将其移动到C盘上，B盘为可以利用的柱子。1.3.3 实现程序function hanta(posfrom,posmiddle,posend,numplate)global nmove%移动次数，调用nmove之前声明nmove为全局变量，且赋值为0if numplate=1 sprintf(从%s移到%s,posfrom,posend) nmove=nmove+1; returnendtry hanta(posfrom,posend,posmiddle,numplate-1) sprintf(从%s移到%s,posfrom,posend) nmove=nmove+1; hanta(posmiddle,posfrom,posend,numplate-1)catch catch errorend实例：当有3个盘子时，该递归程序的输出为：ans =从A移到Cans =从A移到Bans =从C移到Bans =从A移到Cans =从B移到Aans =从B移到Cans =从A移到Cnmove =71.4 案例：商人安全过河问题三名商人各带一个随从乘船渡河，一只小船只能容纳二人，由他们自己划行。随从们密约，在河的任一岸，一旦随从的人数比商人多，就杀人越货。但是如何乘船渡河的大权掌握在商人们手中。商人们怎样才能安全渡河呢？1.4.1 问题分析设商人某岸上的人数为x，随从人数为y；此岸彼岸商人与随从的人数必须满足以下下面其中一种情况：（1） x不等于0时，x=y；（2） x等于0以上两种称为安全状态。每次划行时做出决策：可能的决策必须满足两个条件：（1） 船上载人不超过2人；（2） 确保出发点岸上和到达对岸后，两岸的人数必须同时满足安全状态1.4.2 模型建立设第k次渡河前此岸的商人数为xk，随从人数为yk，k=1,2, xk，yk=0,1,2,3.将二维向量Sk=(xk，yk)定义为状态（即第k-1次渡河后的此岸状态）。安全渡河条件下的状态几何称为允许状态集合，记作S，根据问题分析很容易得到S(x,y)|(0,0),(0,1),(0,2),(0,3),(1,1),(2,2)(3,0),(3,1),(3,2),(3,3) （1）注：为什么不能取(1,0),(2,0),(2,1)?记第k次渡船上的商人数为uk，随从数为vk，将二维向量dk= uk，vk定义为决策。允许决策集合记作D，由小船的容量可知D(u,v)| (1,0),(0,1),(1,1),(2,0),(0,2) （2）因为k为奇数时船从此岸驶向彼岸，k为偶数时船由彼岸驶回此岸，所以状态Sk，随决策dk变化的规律是Sk+1= Sk+(-1)k dk （3）如：（3）式称为状态转移律。这样，制定安全渡河方案归结为如下的多步决策问题：求决策dkD（k=1,2,n），使状态SkS按照转移规律（3）由初始状态S1=(3,3)经过有限步n到达状态Sn+1=(0,0)以上建立了一个简单的模型！但还没有求解。1.4.3 模型求解如何求解上述模型？手工推算？还是借助计算机，通过设计求解算法来求解该问题？1.4.4 进一步的思考题1、 最少应划行多少次船才能过河？2、 总共有多少种过河方法？3、 还有无其他建模方法（还能否建立其他数学模型？比如图论模型）1.4.5 程序运行结果程序名：my_guohe_main（见后面）下面为程序中的注释，通过这个可以理解下面的运行结果输出。思考题：请根据结果输出，将运行结果转化为文字叙述。%用于回溯的矩阵%mat_back = 阶段，状态（决策前），决策，阶段（决策后），状态（决策后），是否结束%mat_back对应列序号： 1 2 3 4 5 6 7 8 9运行结果：ans = 129 9ans = 11 0 2 0 2 0 0 0 1 11 1 1 1 1 0 0 0 1 13 0 2 0 2 0 0 0 1 13 1 1 1 1 0 0 0 1 15 0 2 0 2 0 0 0 1第 1个策略cur_dec_seq = 1 3 3 1 1 2 2 2 0 2 2 2 1 0 3 3 2 0 3 3 2 0 2 4 3 0 0 4 3 0 0 1 5 3 1 0 5 3 1 2 0 6 1 1 0 6 1 1 1 1 7 2 2 0 7 2 2 2 0 8 0 2 0 8 0 2 0 1 9 0 3 0 9 0 3 0 2 10 0 1 0 10 0 1 0 1 11 0 2 0 11 0 2 0 2 0 0 0 1注：这里省略另外4个策略。思考题：通过程序及其注释，分析结果输出的5种过河策略，并翻译为通俗易懂的过河策略。注：对于一组确定的策略（比如上面的5种），实际上“第i行的第7,8列元素应该分别等于第i+1行的第2，3列元素”。1.4.6 递归算法求解程序function my_guohe_main%2004-12-9 调试完毕，共有5种不同的决策clear allglobal fid numans mat_back%用于回溯的矩阵%mat_back = 阶段，状态（决策前），决策，阶段（决策后），状态（决策后），是否结束%mat_back对应列序号： 1 2 3 4 5 6 7 8 9%状态和决策各2个元素存储mat_back=zeros(1,9);%init%是否结束：用1表示结束；0表示未结束；numans = 0 ;%决策序列个数fid = fopen(debug.txt,w);if fid 3 disp(error open file) returnend%调用递归求解程序getd(3 3,1,-1 -1 -1 -1,-1 -1 -1 -1)fprintf(fid,nTotal ans num=%d,numans)fclose(fid)%mat_backsize(mat_back)%回溯得到过河策略mat_back = unique(mat_back,rows);%消除重复行row_no=find(mat_back(:,9)=1);%最后一步决策mat_back(row_no,:)for i=1:length(row_no),%逐个回溯最优解 disp(sprintf(第%4d个策略,i) cur_step_status=mat_back(row_no(i),1:3); step = mat_back(row_no(i),1); %决策序列 cur_dec_seq(step,:)= mat_back(row_no(i),1:9); %只要还没有到初态【3，3】，则继续回溯 while vec_equal(cur_step_status(2:3),3 3), id_mat= find(mat_back(:,9)=1);%?2004-12-9 mid_mat_back= mat_back(id_mat,6 7 8); %oldversion:r,id=vec_in_matrix(cur_step_status,mat_back(:,6 7 8); r,id=vec_in_matrix(cur_step_status,mid_mat_back); id = id_mat(id); %disp(step - 1) step=step-1; if r=1 | id =0 disp(err huisuo) pause return end %决策序列 %size_mat_back=size(mat_back) cur_dec_seq(step,:)=mat_back(id,1:9); %当前 s t e p状态 cur_step_status = mat_back(id,1:3); end%while if step=1 disp(step = 1) cur_step_status cur_dec_seq return end cur_dec_seqend%forfunction getd(status,step,ever1,ever2)%reference：建模示例之二：商人们怎样安全过河%数学模型姜启源%编写一个递归程序，求解此模型%两个重要的变量（可以避免不不要搜索及死循环）%ever1 列分布，第1、2列 状态，第3、4列 决策（此岸已经做的状态、决策）%ever2 列分布，第1、2列 状态，第3、4列 决策（彼岸已经做的状态、决策）% getd(3 3,1,-1 -1 -1 -1,-1 -1 -1 -1,0 0)%status 在此岸人数状态（商人数，随从数）global fid global numans mat_back%if status(1)=0 & status(2)=0 % if status(1:2)=0 0 %成功过河 sprintf(First:Game Over!Sucess! congratulation!) returnend%允许状态allow_s=0 1 %不能包含0 0行，因为前面已经处理了返回 0 2 0 3 3 0 3 1 3 2 3 3 1 1 2 2;%过河允许决策allow_d= 1 0 0 1 1 1 2 0 0 2;row_d,col_d=size(allow_d);%计算可行决策数目，从而一开始就有row_d种决策for i=1:row_d,%for99 cur_status = status + getid(step)*allow_d(i,:);%当前此岸状态 %当前决策是否合法 if cur_status(1) = 3 & cur_status(2) =0 & cur_status(2) =0 & . ( cur_status(1)=0 |cur_status(2)=0| . ( cur_status(1)0 & cur_status(1)=cur_status(2) ) &. (cur_status(1)=3 & cur_status(2)=3)%middle cursd=cur_status,allow_d(i,:);% status decision if cur_status(1)=0 & cur_status(2)=0 %88 很重要的递归结束条件 %mat_back=阶段，状态（决策前），决策，阶段（决策后），状态（决策后），是否结束 mat_back(getrows(mat_back)+1,:)=step,status,allow_d(i,:),0,0,0,1;% numans = numans + 1; sprintf(success over) return else%88 if getid(step)=1 , %77 %if11 if vec_in_matrix(cur_status,allow_s)=1 & vec_in_matrix(cursd,ever1)=0 ever1=ever1;cursd; %mat_back=阶段，状态（决策前），决策，阶段（决策后），状态（决策后），是否结束 mat_back(getrows(mat_back)+1,:)=step,status,allow_d(i,:),step+1,cur_status,0; getd(cur_status,step+1,ever1,ever2); end%if11 else %if22 if vec_in_matrix(cur_status,allow_s)=1 & vec_in_matrix(cursd,ever2)=0 ever2=ever2;cursd; %mat_back=阶段，状态（决策前），决策，阶段（决策后），状态（决策后），是否结束 mat_back(getrows(mat_back) +1,:)=step,status,allow_d(i,:),step+1,cur_status,0; getd(cur_status,step+1,ever1,ever2); end%if 22 end %if 77 end%if 88 else end %if middle end%for 99function r,id=vec_in_matrix(cur,matrix)%判定向量是否在矩阵中row,col=size(matrix);for i=1:row if vec_equal(cur,matrix(i,:)=1 r=1; id=i;%行号 return end endr=0;function r=vec_equal(vec1,vec2)r=0;%判定向量是否相等if vec1=vec2, r=1;endfunction r=getid(step)%根据当前决策步骤step的奇偶性判定if mod(step,2)=1%（奇数）此岸-彼岸 r=-1;else r=1;%(偶数)此岸-彼岸endfunction rows= getrows(mat)rows,tmp_tmp=size(mat);%返回函数2 优化模型及其求解2.1 案例：背包问题有一组物品S，共有9件，其中第i件重 ，价值，从S中取出一些物品出来装背包，使总价值最大，而不超过总重量的给定上限30kg。 i123456789(kg)2112.5106543(元)104530100150902001803002.1.1 问题分析这是一个典型的最优化问题，优化目标是总价值最大，决策是决定装哪些物品，而装载物品又受到背包所能承受重量30kg的限制。因此可以建立该问题的最优化数学模型，而且是01整数规划模型。2.1.2 变量与符号说明：用来表示是否装载第i件物品，如果表示不装载该物品，如果表示装载该物品（），令。：第i件物品重量（），单位：kg。：第i件物品价值（），单位：元。：表示按照决策，得到的装载物品的总价值，单位：元。2.1.3 模型建立通过前面的分析及变量的定义，可得。原问题抽象为如下的01整数线性规划模型：2.1.4 模型求解及结果可以采用求解线性规划的Lindo软件求解，下面采用求解可以求解非线性规划的Lingo软件求解。Lingo功能非常强大，其语法本身较复杂，这里只用到最简单的语法。直接参看代码。Lindo基本语法（1） 模型用“MODEL:”开始，“END”结束，参考下面的Lingo程序；（2） 语句用分号中断。注释行用感叹号开头，通样是分号结束。（3） 如果变量只取0或1，则用BIN设置，比如BIN(x2)，表示变量x2只能取0或1。（4） 如果变量只取整数，则用GIN设置，比如GIN(x2)，表示变量x2为整数变量。Lingo程序如下：!求解背包问题的线性整数规划模型（01整数规划）;MODEL:!目标函数; max = 10*x1+ 45*x2+30*x3+100*x4 +150*x5+90*x6+200*x7+180*x8+300*x9;!重量约束;2*x1+x2+x3+2.5*x4 +10*x5 +6*x6+ 5*x7 + 4*x8+ 3*x9=30;!整数约束的设置;BIN(x1);BIN(x2);BIN(x3);BIN(x4);BIN(x5);BIN(x6);BIN(x7);BIN(x8);BIN(x9);END求解显示结果：Global optimal solution found at iteration: 0 Objective value: 1015.000 Variable Value Reduced Cost X1 1.000000 -10.00000 X2 1.000000 -45.00000 X3 1.000000 -30.00000 X4 1.000000 -100.0000 X5 1.000000 -150.0000 X6 0.000000 -90.00000 X7 1.000000 -200.0000 X8 1.000000 -180.0000 X9 1.000000 -300.0000 Row Slack or Surplus Dual Price 1 1015.000 1.000000 2 1.500000 0.000000结果显示已经找到了最优解。通过结果可见：除了不装载第6件物品外，其余全部装载，这种装包方法达到总价值最大为1015元。下面介绍其它求解方法。2.1.5 贪婪法建立优选指标，即第i件物品单位重量的价值（元/kg）。选择装背包的物品时，就从优选指标中最大的开始选，直到不能再装为止。由已知数据，可得i1234567895453040151540451002.1.6 贪婪法求解程序%背包问题求解得贪婪算法实现w=2112.5106543;v=10453010015090200180300;g=v./w; %单位质量的价值（指标）maxw=30;%对指标排序，以便根据指标选择物品y,idx=sort(g);curw =0;%存储当前装的总重curv =0;%存储当前装的价值hasidx=;%存储当前装的物品序号for i=length(g):-1:1, %从指标g最大的开始选择 id = idx(i);%important if curw + w(id) = maxw, curw = curw + w(id); curv = curv + v(id); hasidx = hasidx ,id; else break;%结束，不能再装 endenddisp(sprintf(总共装载物品数=%5d件,length(hasidx)disp(sprintf(质量=%10.2fkg,curw)disp(sprintf(价值=%10.2f元,sum(v(hasidx)%实际上curw = sum(w(hasidx)% curv =sum(v(hasidx)2.1.7 贪婪法求解结果运行程序输出为：总共装载物品数= 7件质量= 22.50kg价值= 945.00元按照贪婪法策略，得到装载物品的顺序依次为：9，8，2，7，4，3，6。思考：以上装载结果是否达到最优？2.1.8 穷举法求解程序%背包问题求解得穷举法得到最优解的程序w=2112.5106543;%各件物品重量v=10453010015090200180300;%价值maxw=30;%背包承受重量上限optvalue = -1; %最优目标（最大价值）初始化for x1=0:1, for x2=0:1, for x3=0:1, for x4=0:1, for x5=0:1, for x6=0:1, for x7=0:1, for x8=0:1, for x9=0:1, x=x1 x2 x3 x4 x5 x6 x7 x8 x9;%决策向量 %不超过重量限制，也要比以前装包方案好 if (dot(x,w)optvalue), optx=x;%存储 optvalue=dot(x,v);%存储 end end end end end end end end endendhasidx = find(optx=0);optxdisp(sprintf(总共装载物品数=%5d件,length(hasidx)disp(sprintf(质量=%10.2fkg,dot(optx,w)disp(sprintf(价值=%10.2f元,dot(optx,v) %实际上：dot(optx,v)=sum(v(hasidx)%dot(optx,w)=sum(w(hasidx)2.1.9 穷举法程序运行结果总共装载物品数= 8件质量= 28.50kg价值= 1015.00元hasidx =1 2 3 4 5 7 8 9只有编号为6的物品没有装上。2.2 案例：高速公路问题A城和B城之间准备建一条高速公路，B城位于A城正南20公里和正东30公里交汇处，它们之间有东西走向连绵起伏的山脉。公路造价与地形特点有关，图4.2.4给出了整个地区的大致地貌情况，显示可分为三条沿东西方向的地形带。你的任务是建立一个数学模型，在给定三种地形上每公里的建造费用的情况下，确定最便宜的路线。图中直线AB显然是路径最短的，但不一定最便宜。而路径ARSB过山地的路段最短，但是否是最好的路径呢？你怎样使你的模型适合于下面两个限制条件的情况呢？1. 当道路转弯是，角度至少为1400。2. 道路必须通过一个已知地点（如P）。 平原 R P 高地高山高地S 平原A B2.2.1 问题分析在建设高速公路时，总是希望建造费用最小。如果要建造的起点、终点在同一地貌中，那么最佳路线则是两点间连接的线段，这样费用则最省。因此本问题是一个典型的最优化问题，以建造费用最小为目标，需要作出的决策则是确定在各个地貌交界处的汇合点。2.2.2 变量说明：在 第i个汇合点上的横坐标（以左下角为直角坐标原点），I1，2，。，4；x530（指目的地B点的横坐标），设；：第i段南北方向的长度（i1，2，。，5）；Si：在第i段上地所建公路的长度（i1，2，。，5）；由问题分析可知，C1：平原每公里的造价（单位：万元/公里）C2：高地每公里的造价（单位：万元/公里）C3：高山每公里的造价（单位：万元/公里）2.2.3 模型假设1、 假设在相同地貌中修改高速公路，建造费用与公路长度成正比；2、 假设在相同地貌中修改高速为直线。在理论上，可以使得建造费用最少，当然实际中一般达不到。2.2.4 模型建立在A城与B城之间建造一条高速公路的问题可以转化为下面的非线性规划模型。优化目标是在A城与B城之间建造高速公路的费用。2.2.5 模型求解这里采用Matlab编程求解。模型求解时，分别去Ci如下。平原每公里的造价C1：400（单位：万元）高地每公里的造价C2：800（单位：万元）高山每公里的造价C3：1200（单位：万元）（注意：实际建模时必须查找资料来确定参数或者题目给定有数据）输入主程序model_p97.m，运行结果如下：model_p97optans = 2.2584e+004len = 38.9350ans = 12.1731 14.3323 15.6677 17.8269求解程序见附录。2.2.6 模型结果及分析通过求解可知，为了使得建造费用最小。建造地点的选择宜采取下列结果。建造总费用为2.2584亿元。总长度为38.9350公里。2.2.7 求解模型的主程序文件程序名：model_p97function x=model_p97%数学建模教材 P97 高速公路clear allglobal C LC=400 800 1200;L=4 4 4 4 4;x=fmincon(objfun_97,1,1,1,1,zeros(1,4),ones(1,4)*30,mycon_p97);optans=objfun_97(x)C=ones(3,1);len = objfun_97(x)2、模型中描述目标函数的Matlab程序objfun_97.mfunction obj=objfun_97(x)global C Lobj= C(1)*sqrt(L(1)2+x(1)2) + C(2)*sqrt(L(2)2+(x(2)-x(1)2) + . C(3)*sqrt(L(3)2+(x(3)-x(2)2) + C(2)*sqrt(L(4)2+(x(4)-x(3)2)+. C(1)*sqrt(L(5)2+(30-x(4)2);3、模型中描述约束条件的Matlab函数mycon_p97.mfunction c,ceq=mycon_p97(x)%c(1)=x(1)-x(2);c(2)=x(2)-x(3);c(3)=x(3)-x(4);c(4)=x(4)-30;ceq=;2.3 随机跳跃法2.3.1 随机跳跃法简介对于如下的优化模型，可以采用随机跳跃法求解。思路：产生若干组0,1均匀分布的随机数，然后求出对应每组随机数的X，随机产生m组这样的点X，然后从这些点中找出目标函数的最小值，则找到一个近似最优解。特点：1、 方法简单，但并不适用于变量很多的情形。如果不考虑运行效率，倒是可以用该方法一试。随机产生的点越多，就越接近于最优解。2、 一般不用于实时系统。2.3.2 求解高速公路问题的随机跳跃法程序程序：solvP97rand.m%2004-12-3%随机跳跃法示例：求解高速公路问题%global C LC=400 800 1200;L=4 4 4 4 4;time_begin = clock;N=4;%决策变量个数 TotalCount=1e5; a=0; b=30; opt_obj = inf;%init opt_dec = zeros(1,4); for c=1:TotalCount, for i=1:N, x(i)=a + rand*(b-a); end %for i cur_obj = objfun_97(x); if cur_obj opt_obj, opt_obj = cur_obj; opt_dec = x; end end %for c disp(sprintf(最优目标%20.6f,opt_obj) for i=1:N, disp(sprintf(x(%d)=%20.4f,i,x(i) end exp_time =etime(clock,time_begin)C=ones(3,1);%相对于讲目标函数中单价设为1，则计算结果应为总长度S1+S2+S3+S4+S5len = objfun_97(opt_dec)2.3.3 程序运行结果运用随机跳跃法搜索高速问题的近似最优解，程序输出结果如下：最优目标 22623.036913x(1)= 18.1229x(2)= 0.5607x(3)= 28.1015x(4)= 24.4806exp_time = 15.5500len = 38.6446发现超出最优值，可见已经非常接近最优值了，这里采用随机跳跃法的效果较好。思考：为什么随机跳跃法在这里运用的效果这么好？2.4 网格法2.4.1 网格法简介对于如下的优化模型，也可以采用网格法求解。思路：将各决策变量的取值区间等分为多个区间，比如m等分，而在搜索中对各决策变量的取值，只从这等分点上取出组合，各变量有m+1个可能的取值。由于m等分后，各决策变量就只搜索m+1个点，模型有n个决策变量，则要搜索个可行解。然后从这些点中找出目标函数的最小值，则找到一个近似最优解。特点：1、 方法简单，并不适用于变量很多的情形。如果不考虑运行效率，倒是可以用该方法一试。2、 如果要结果更精确，m可以取更大些，就越接近于最优解，当如搜索点更多，更费时。3、 一般不用于实时系统。2.4.2 求解高速公路问题的网格法程序%2004-12-3%网格法示例：求解高速公路问题%clear allglobal C LC=400 800 1200;L=4 4 4 4 4;time_begin = clock;N=4;%决策变量个数a=0; b=30; opt_obj = inf;%init opt_dec = zeros(1,4);gridnum=30; step=b/gridnum count=0; for x1=0:step: b, for x2=0:step:b, for x3=0:step:b, for x4=0:step:b, x = x1 x2 x3 x4; count=count+1; cu