下载文档原格式
重庆大学（商仆过河模型）开课课程：数学模型指导教师：黄光辉小组成员：庄楚斌（20134760）自动化07班张俊铭（20133227）材料加工01班吴慧（20131966）数学01班时间：2015年3月8日一、问题提出3名商人带3名随从乘一条小船过河，小船每次只能承载至多两人。
随从们密约, 在河的任一岸, 一旦随从的人数比商人多,就杀人越货.乘船渡河的方案由商人决定，商人们如何才能安全渡河呢？二、问题分析商人与随从过河问题可以视为是一个多步决策的过程，通过多次优化，从而得到一个全局最优的决策方案。
决策的每一步，即船从此岸到达彼岸，都要对船上的商人和仆人数做出决策。
在保证河的任一岸均有商人数比随从人数多和小船每次最多只能承载两人的前提下，经有限步使所有人员到达彼岸。
三、模型假设商人和随从都会划船，天气很好，无大风大浪，且船的质量很好，可以保证很多次安全的运载商人和随从。
四、模型建立x~第k次渡河前此岸的商人数，k y~第k次渡河前此岸的随从数kx, k y=0,1,2,3; k=1,2,……kS=(k x, k y, c k )~过程的状态，其中k x, k y, c k 分别表示对应时刻此k岸的商人，仆人数以及船的行进方向，其中c 取值1表示即将向彼岸运行，为0表示即将向此岸运行S ~ 允许状态集合，S={(x , y )| x =0, y =0,1,2,3; x =3 ,y =0,1,2,3; x =y =1,2} k u ~第k 次渡船上的商人数k v ~第k 次渡船上的随从数k d =(k u , k v )~决策，D={(u , v )| 21≤+≤v u ,k u , k v =0,1,2} ~允许决策集合k =1,2,… …因为k 为奇数时船从此岸驶向彼岸，k 为偶数时船从彼岸驶向此岸，所以状态k S 随决策k d 的变化规律是1+k S =k S +k )1(-k d ~状态转移律求k d ∈D(k =1,2, …n), 使k S ∈S, 并按转移律由1S =(3,3,1)到达状态1+n S =(0,0,0(1))。
function [pTS,fmin]=grTravSale(C)% Function [pTS,fmin]=grTravSale(C) solve the nonsymmetrical% traveling salesman problem.% Input parameter:% C(n,n) - matrix of distances between cities,% maybe, nonsymmetrical;% n - number of cities.% Output parameters:% pTS(n) - the order of cities;% fmin - length of way.% Uses the reduction to integer LP-problem:% Look: Miller C.E., Tucker A. W., Zemlin R. A.% Integer Programming Formulation of Traveling Salesman Problems. % J.ACM, 1960, Vol.7, p. 326-329.% Needed other products: MIQP.M.% This software may be free downloaded from site:% http://control.ee.ethz.ch/~hybrid/miqp/% Author: Sergiy Iglin% e-mail: siglin@yandex.ru% personal page: http://iglin.exponenta.ru% ============= Input data validation ==================if nargin<1,error('There are no input data!')endif ~isnumeric(C),error('The array C must be numeric!')endif ~isreal(C),error('The array C must be real!')ends=size(C); % size of array Cif length(s)~=2,error('The array C must be 2D!')endif s(1)~=s(2),error('Matrix C must be square!')endif s(1)<3,error('Must be not less than 3 cities!')end% ============ Size of problem ====================n=s(1); % number of vertexesm=n*(n-1); % number of arrows% ============ Parameters of integer LP problem ========Aeq=[]; % for the matrix of the boundary equationsfor k1=1:n,z1=zeros(n);z1(k1,:)=1;z2=[z1;eye(n)];Aeq=[Aeq z2([1:2*n-1],setdiff([1:n],k1))];endAeq=[Aeq zeros(2*n-1,n-1)];A=[]; % for the matrix of the boundary inequationsfor k1=2:n,z1=[];for k2=1:n,z2=eye(n)*(n-1)*(k2==k1);z1=[z1 z2(setdiff([2:n],k1),setdiff([1:n],k2))];endz2=-eye(n);z2(:,k1)=z2(:,k1)+1;A=[A;[z1 z2(setdiff([2:n],k1),2:n)]];endbeq=ones(2*n-1,1); % the right parts of the boundary equations b=ones((n-1)*(n-2),1)*(n-2); % the right parts of the boundary inequationsC1=C'+diag(ones(1,n)*NaN);C2=C1(:);c=[C2(~isnan(C2));zeros(n-1,1)]; % the factors for objective functionvlb=[zeros(m,1);-inf*ones(n-1,1)]; % the lower boundsvub=[ones(m,1);inf*ones(n-1,1)]; % the upper boundsH=zeros(n^2-1); % Hessian% ============= We solve the MIQP problem ========== [xmin,fmin]=MIQP(H,c,A,b,Aeq,beq,[1:m],vlb,vub);% ============= We return the results ==============eik=round(xmin(1:m)); % the arrows of the waye1=[zeros(1,n-1);reshape(eik,n,n-1)];e2=[e1(:);0]; % we add zero to a diagonale3=(reshape(e2,n,n))'; % the matrix of the waypTS=[1 find(e3(1,:))]; % we build the waywhile pTS(end)>1, % we add the city to the waypTS=[pTS find(e3(pTS(end),:))];endreturn。
走遍全中国方案的研究摘要本文通过对走遍中国各省会城市、直辖市和港澳台的最优路径选择问题进行分析，发现这是一个十分典型的旅行商问题（Traveling Salseman Problem ），即寻找一条遍历n 个城市（在本文中为34个城市）的最短路径。
我们搜索了这34个城市的经纬度和部分列车、航班时刻表等各方面信息，综合省钱、省时、方便等因素，进一步深入并细化，从而得到判断各订票方案的准则。
针对问题一，我们利用欧几里得平面知识，由公式0002901800290A B x R A y R ππ+⎧=⋅⎪⎪⎨-⎪=⎪⎩将该34个城市的经纬度转化为平面直角坐标，从而将其简化成二维平面问题。
目前，解旅行商问题（tsp 问题）的方法有许多种，本文采用了较为先进的遗传算法。
遗传算法是目前解决组合优化问题最有效的工具之一，本文介绍了遗传算法的基本原理，讨论了遗传算法中的有关遗传算子设计等方面的技术。
由于该算法在搜索空间中同时考虑了许多点，这样就减少了收敛于局部极小的可能，也增加了处理的并行性。
同时，我们利用MATLAB 软件编辑了相关程序，计算出了遍历该34个城市的最短路径，其路径长度为14661km 。
对于问题二，在只考虑旅行路线最经济的前提下，结合第一问得出的最优路径，我们收集了这34个城市的列车和航班时刻表等信息，从而找出最经济的订票方案，并得其花费为11426元。
针对问题三，在综合考虑省时、省钱、方便等诸多因素的前提下制定订票方案的评价准则，我们运用了层次分析法对其进行研究。
根据对旅行过程中省时，省钱及方便的偏重程度，我们相应地给出了判断矩阵111571513731⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦A ，然后对其进行一致性检验，发现其不一致程度在容许范围内。
因此我们利用其最大特征根max λ对应的特征向量w 作为比较因素的权向量，并得到以下评价表达式：12120.0740.2830.643*(0.8330.167)C x x x x =+++。
1.问题基本描述:求一个旅行商经过N个城市最后回到出发点的最短路径.即,在一个无向带权图的邻接矩阵中,求一个最短环包括所有顶点.2.解法:1)动态规划:假设从顶点i出发，令d(i,V’)表示从顶点i出发经过V’中各个顶点一次且仅一次，最后回到出发点i的最短路径的长度，开始时，V’=V-{i}，于是，旅行商问题的动态规划函数为：d(i,V’) = min{c ik + d(k,V’-{k})} (k∈V’)1)d(k,{}) = c ki (k ≠ i)2)简单来说,就是用递归表达:从出发点0到1号点,假设1是第一个,则剩下的路程就是从1经过剩下的点最后回到0点的最短路径. 所以当V’为空的时候, d(k,{}) = c ki (k ≠ i), 找的是最后一个点到0点的距离.递归求解1之后,再继续求V’之中剩下的点,最后找出min.如果按照这个思想直接做,对于每一个i都要递归剩下的V中所有的点,所以这样的时间复杂度就近似于N!,其中有很多重复的工作.可以从小的集合到大的集合算,并存入一个二维数组,这样当加入一个节点时,就可以用到之前的结果,如四个点的情况:动态填表:表中元素代表第i个节点经过V集合中的点最后到0点的最短值.如果有多个值,取其中最小的一个.这样一共循环(2^(N-1)-1)*(N-1)次,就把时间复杂度缩小到O(N*2)的级别.核心伪代码如下:{for (i =1;i<n;i++) //初始化第0列d[i][0]=c[i][0];for( j=1;j<2^(N-1)-1;j++)for(i=1 ; i<n ;i++){if(子集Vj中不包含i){对Vj中的每个元素k，计算d[i][Vj] = min{c[i][k] + d[k][{Vj-k}] | 每一个k∈Vj};}}对V[2^(n-1)-1]中的每个元素k,计算:d[0][2^(n-1)-1] = min{c[0][k] + d[k][2^(n-1)-2]};输出最短路径:d[0][2^(n-1)-1];}。
function youhuafunD=code;N=50; % Tunablemaxgen=50; % Tunablecrossrate=0.5; %Tunablemuterate=0.08; %Tunablegeneration=1;num = length(D);fatherrand=randint(num,N,3);score = zeros(maxgen,N);while generation<=maxgenind=randperm(N-2)+2; % 随机配对交叉A=fatherrand(:,ind(1:(N-2)/2));B=fatherrand(:,ind((N-2)/2+1:end));% 多点交叉rnd=rand(num,(N-2)/2);ind=rnd tmp=A(ind);A(ind)=B(ind);B(ind)=tmp;% % 两点交叉% for kk=1:(N-2)/2% rndtmp=randint(1,1,num)+1;% tmp=A(1:rndtmp,kk);% A(1:rndtmp,kk)=B(1:rndtmp,kk);% B(1:rndtmp,kk)=tmp;% endfatherrand=[fatherrand(:,1:2),A,B];% 变异rnd=rand(num,N);ind=rnd [m,n]=size(ind);tmp=randint(m,n,2)+1;tmp(:,1:2)=0;fatherrand=tmp+fatherrand;fatherrand=mod(fatherrand,3);% fatherrand(ind)=tmp;%评价、选择scoreN=scorefun(fatherrand,D);% 求得N个个体的评价函数score(generation,:)=scoreN;[scoreSort,scoreind]=sort(scoreN);sumscore=cumsum(scoreSort);sumscore=sumscore./sumscore(end);childind(1:2)=scoreind(end-1:end);for k=3:Ntmprnd=rand;tmpind=tmprnd difind=[0,diff(t mpind)];if ~any(difind)difind(1)=1;endchildind(k)=scoreind(logical(difind));endfatherrand=fatherrand(:,childind);generation=generation+1;end% scoremaxV=max(score,[],2);minV=11*300-maxV;plot(minV,'*');title('各代的目标函数值');F4=D(:,4);FF4=F4-fatherrand(:,1);FF4=max(FF4,1);D(:,5)=FF4;save DData Dfunction D=codeload youhua.mat% properties F2 and F3F1=A(:,1);F2=A(:,2);F3=A(:,3);if (max(F2)>1450)||(min(F2)<=900)error('DATA property F2 exceed it''s range(900,1450]')end% get group property F1 of data, according to F2 value F4=zeros(size(F1));for ite=11:-1:1index=find(F2<=900+ite*50);F4(index)=ite;endD=[F1,F2,F3,F4];function ScoreN=scorefun(fatherrand,D)F3=D(:,3);F4=D(:,4);N=size(fatherrand,2);FF4=F4*ones(1,N);FF4rnd=FF4-fatherrand;FF4rnd=max(FF4rnd,1);ScoreN=ones(1,N)*300*11;% 这里有待优化for k=1:NFF4k=FF4rnd(:,k);for ite=1:11F0index=find(FF4k==ite);if ~isempty(F0index)tmpMat=F3(F0index);tmpSco=sum(tmpMat);ScoreBin(ite)=mod(tmpSco,300);endendScorek(k)=sum(ScoreBin);endScoreN=ScoreN-Scorek;遗传算法实例:% 下面举例说明遗传算法 %% 求下列函数的最大值 %% f(x)=10*sin(5x)+7*cos(4x) x∈[0,10] %% 将 x 的值用一个10位的二值形式表示为二值问题，一个10位的二值数提供的分辨率是每为 (10-0)/(2^10-1)≈0.01 。
MATLAB多旅行商问题源代码function varargout = mtspf_ga(xy,dmat,salesmen,min_tour,pop_size,num_iter,show_prog,show_res)% MTSPF_GA Fixed Multiple Traveling Salesmen Problem (M-TSP) Genetic Algorithm (GA)% Finds a (near) optimal solution to a variation of the M-TSP by setting% up a GA to search for the shortest route (least distance needed for% each salesman to travel from the start location to individual cities% and back to the original starting place)%% Summary:% 1. Each salesman starts at the first point, and ends at the first% point, but travels to a unique set of cities in between% 2. Except for the first, each city is visited by exactly one salesman%% Note: The Fixed Start/End location is taken to be the first XY point%% Input:% XY (float) is an Nx2 matrix of city locations, where N is the number of cities% DMAT (float) is an NxN matrix of city-to-city distances or costs% SALESMEN (scalar integer) is the number of salesmen to visit the cities% MIN_TOUR (scalar integer) is the minimum tour length for any of the% salesmen, NOT including the start/end point% POP_SIZE (scalar integer) is the size of the population (should be divisible by 8)% NUM_ITER (scalar integer) is the number of desired iterations for the algorithm to run% SHOW_PROG (scalar logical) shows the GA progress if true% SHOW_RES (scalar logical) shows the GA results if true%% Output:% OPT_RTE (integer array) is the best route found by the algorithm% OPT_BRK (integer array) is the list of route break points (these specify the indices% into the route used to obtain the individual salesman routes)% MIN_DIST (scalar float) is the total distance traveled by the salesmen%% Route/Breakpoint Details:% If there are 10 cities and 3 salesmen, a possible route/break% combination might be: rte = [5 6 9 4 2 8 10 3 7], brks = [3 7]% Taken together, these represent the solution [1 5 6 9 1][1 4 2 8 1][1 10 3 7 1],% which designates the routes for the 3 salesmen as follows:% . Salesman 1 travels from city 1 to 5 to 6 to 9 and back to 1% . Salesman 2 travels from city 1 to 4 to 2 to 8 and back to 1% . Salesman 3 travels from city 1 to 10 to 3 to 7 and back to 1%% 2D Example:% n = 35;% xy = 10*rand(n,2);% salesmen = 5;% min_tour = 3;% pop_size = 80;% num_iter = 5e3;% a = meshgrid(1:n);% dmat = reshape(sqrt(sum((xy(a,:)-xy(a',:)).^2,2)),n,n);% [opt_rte,opt_brk,min_dist] = mtspf_ga(xy,dmat,salesmen,min_tour, ... % pop_size,num_iter,1,1);%% 3D Example:% n = 35;% xyz = 10*rand(n,3);% salesmen = 5;% min_tour = 3;% pop_size = 80;% num_iter = 5e3;% a = meshgrid(1:n);% dmat = reshape(sqrt(sum((xyz(a,:)-xyz(a',:)).^2,2)),n,n);% [opt_rte,opt_brk,min_dist] = mtspf_ga(xyz,dmat,salesmen,min_tour, ... % pop_size,num_iter,1,1);%% See also: mtsp_ga, mtspo_ga, mtspof_ga, mtspofs_ga, mtspv_ga, distmat%% Author: Joseph Kirk%Email:*******************% Release: 1.3% Release Date: 6/2/09% Process Inputs and Initialize Defaultsnargs = 8;for k = nargin:nargs-1switch kcase 0xy = 10*rand(40,2);case 1N = size(xy,1);a = meshgrid(1:N);dmat = reshape(sqrt(sum((xy(a,:)-xy(a',:)).^2,2)),N,N);case 2salesmen = 5;case 3min_tour = 2;case 4pop_size = 80;case 5num_iter = 5e3;case 6show_prog = 1;case 7show_res = 1;otherwiseendend% Verify Inputs[N,dims] = size(xy);[nr,nc] = size(dmat);if N ~= nr || N ~= ncerror('Invalid XY or DMAT inputs!')endn = N - 1; % Separate Start/End City% Sanity Checkssalesmen = max(1,min(n,round(real(salesmen(1)))));min_tour = max(1,min(floor(n/salesmen),round(real(min_tour(1))))); pop_size = max(8,8*ceil(pop_size(1)/8));num_iter = max(1,round(real(num_iter(1))));show_prog = logical(show_prog(1));show_res = logical(show_res(1));% Initializations for Route Break Point Selectionnum_brks = salesmen-1;dof = n - min_tour*salesmen; % degrees of freedom addto = ones(1,dof+1);for k = 2:num_brksaddto = cumsum(addto);endcum_prob = cumsum(addto)/sum(addto);% Initialize the Populationspop_rte = zeros(pop_size,n); % population of routespop_brk = zeros(pop_size,num_brks); % population of breaksfor k = 1:pop_sizepop_rte(k,:) = randperm(n)+1;pop_brk(k,:) = randbreaks();end% Select the Colors for the Plotted Routesclr = [1 0 0; 0 0 1; 0.67 0 1; 0 1 0; 1 0.5 0];if salesmen > 5clr = hsv(salesmen);end% Run the GAglobal_min = Inf;total_dist = zeros(1,pop_size);dist_history = zeros(1,num_iter);tmp_pop_rte = zeros(8,n);tmp_pop_brk = zeros(8,num_brks);new_pop_rte = zeros(pop_size,n);new_pop_brk = zeros(pop_size,num_brks);if show_progpfig = figure('Name','MTSPF_GA | Current Best Solution','Numbertitle','off'); endfor iter = 1:num_iter% Evaluate Members of the Populationfor p = 1:pop_sized = 0;p_rte = pop_rte(p,:);p_brk = pop_brk(p,:);rng = [[1 p_brk+1];[p_brk n]]';for s = 1:salesmend = d + dmat(1,p_rte(rng(s,1))); % Add Start Distancefor k = rng(s,1):rng(s,2)-1d = d + dmat(p_rte(k),p_rte(k+1));endd = d + dmat(p_rte(rng(s,2)),1); % Add End Distanceendtotal_dist(p) = d;end% Find the Best Route in the Population[min_dist,index] = min(total_dist);dist_history(iter) = min_dist;if min_dist < global_minglobal_min = min_dist;opt_rte = pop_rte(index,:);opt_brk = pop_brk(index,:);rng = [[1 opt_brk+1];[opt_brk n]]';if show_prog% Plot the Best Routefigure(pfig);for s = 1:salesmenrte = [1 opt_rte(rng(s,1):rng(s,2)) 1];if dims == 3, plot3(xy(rte,1),xy(rte,2),xy(rte,3),'.-','Color',clr(s,:));else plot(xy(rte,1),xy(rte,2),'.-','Color',clr(s,:)); endtitle(sprintf('Total Distance = %1.4f, Iteration = %d',min_dist,iter));hold onendif dims == 3, plot3(xy(1,1),xy(1,2),xy(1,3),'ko');else plot(xy(1,1),xy(1,2),'ko'); endhold offendend% Genetic Algorithm Operatorsrand_grouping = randperm(pop_size);for p = 8:8:pop_sizertes = pop_rte(rand_grouping(p-7:p),:);brks = pop_brk(rand_grouping(p-7:p),:);dists = total_dist(rand_grouping(p-7:p));[ignore,idx] = min(dists);best_of_8_rte = rtes(idx,:);best_of_8_brk = brks(idx,:);rte_ins_pts = sort(ceil(n*rand(1,2)));I = rte_ins_pts(1);J = rte_ins_pts(2);for k = 1:8 % Generate New Solutionstmp_pop_rte(k,:) = best_of_8_rte;tmp_pop_brk(k,:) = best_of_8_brk;switch kcase 2 % Fliptmp_pop_rte(k,I:J) = fliplr(tmp_pop_rte(k,I:J));case 3 % Swaptmp_pop_rte(k,[I J]) = tmp_pop_rte(k,[J I]);case 4 % Slidetmp_pop_rte(k,I:J) = tmp_pop_rte(k,[I+1:J I]);case 5 % Modify Breakstmp_pop_brk(k,:) = randbreaks();case 6 % Flip, Modify Breakstmp_pop_rte(k,I:J) = fliplr(tmp_pop_rte(k,I:J));tmp_pop_brk(k,:) = randbreaks();case 7 % Swap, Modify Breakstmp_pop_rte(k,[I J]) = tmp_pop_rte(k,[J I]);tmp_pop_brk(k,:) = randbreaks();case 8 % Slide, Modify Breakstmp_pop_rte(k,I:J) = tmp_pop_rte(k,[I+1:J I]);tmp_pop_brk(k,:) = randbreaks();otherwise % Do Nothingendendnew_pop_rte(p-7:p,:) = tmp_pop_rte;new_pop_brk(p-7:p,:) = tmp_pop_brk;endpop_rte = new_pop_rte;pop_brk = new_pop_brk;endif show_res% Plotsfigure('Name','MTSPF_GA | Results','Numbertitle','off');subplot(2,2,1);if dims == 3, plot3(xy(:,1),xy(:,2),xy(:,3),'k.');else plot(xy(:,1),xy(:,2),'k.'); endtitle('City Locations');subplot(2,2,2);imagesc(dmat([1 opt_rte],[1 opt_rte]));title('Distance Matrix');subplot(2,2,3);rng = [[1 opt_brk+1];[opt_brk n]]';for s = 1:salesmenrte = [1 opt_rte(rng(s,1):rng(s,2)) 1];if dims == 3, plot3(xy(rte,1),xy(rte,2),xy(rte,3),'.-','Color',clr(s,:));else plot(xy(rte,1),xy(rte,2),'.-','Color',clr(s,:)); endtitle(sprintf('Total Distance = %1.4f',min_dist));hold on;endif dims == 3, plot3(xy(1,1),xy(1,2),xy(1,3),'ko');else plot(xy(1,1),xy(1,2),'ko'); endsubplot(2,2,4);plot(dist_history,'b','LineWidth',2);title('Best Solution History');set(gca,'XLim',[0 num_iter+1],'YLim',[0 1.1*max([1 dist_history])]);end% Return Outputsif nargoutvarargout{1} = opt_rte;varargout{2} = opt_brk;varargout{3} = min_dist;end% Generate Random Set of Break Pointsfunction breaks = randbreaks()if min_tour == 1 % No Constraints on Breakstmp_brks = randperm(n-1);breaks = sort(tmp_brks(1:num_brks));else % Force Breaks to be at Least the Minimum Tour Length num_adjust = find(rand < cum_prob,1)-1;spaces = ceil(num_brks*rand(1,num_adjust));adjust = zeros(1,num_brks);for kk = 1:num_brksadjust(kk) = sum(spaces == kk);endbreaks = min_tour*(1:num_brks) + cumsum(adjust);endendend。
