2021年数维杯国际大学生数学建模C题冠状病毒下的旅游业未来发展规划求解全过程文档及程序

2021年数维杯国际大学生数学建模

C题冠状病毒下的旅游业未来发展规划

原题再现：

旅游业是一个具有高度关联性的复合型产业。它不仅与交通运输业、餐饮业、旅游服务业直接相关，而且与第三产业的大部分行业密切相关。旅游业带动的消费支出主要包括三部分：旅游者在旅游目的地购买商品和服务的支出，旅游过程中产生的交通、住宿、餐饮、娱乐、休闲等服务的支出，旅游景区建设和旅游企业各项基础设施建设费用等。
然而，新冠病毒大流行对各行各业都产生了不同程度的影响。这对旅游业有较大的负面影响。一些国家有意或无意地采取措施防止这一流行病的蔓延，严重阻碍了旅游业的发展。旅游业的萧条必然会影响旅游者的劳动欲望，影响旅游企业和政府部门的收入，进而导致经济停滞或回归。在新冠病毒的背景下，一些国家也主动探索新的旅游模式，但认为已经达成了普遍共识。对于我国这样一个人口众多、假期长短不一的国家，如何在疫情背景下权衡各方利益，进而提出科学的旅游发展规划显得非常有意义。
附件1给出了北京部分景点的名称和经纬度信息。请结合数学建模的相关理论，解决以下五个问题：
问题1：参考附件1中的数据。请收集全国各省景区的基础数据。特别需要收集目前景区的评价水平和每天的最大接待量，直观展示这些景区的分布特征。
问题2：能否利用收集到的景点和接待能力，对这些景点进行定量评价和合理分类？
问题3：在疫情背景下，目前各景区的限价已经成为一种新的模式。能否根据不同的流行程度为景区提供合理的限流模型，并以某景区为例给出具体的定量分析结果？
问题4：能否运用数学建模方法，提出同时考虑防疫、旅游收入、旅游体验的定量模型，并对预期效果进行仿真分析。
问题5：能否为政府管理部门提供疫情期间不同景区不超过2页的差异化管理方案？

整体求解过程概述(摘要)

COVID-19对旅游业产生了不同程度的影响。必须采取有效措施解决这些问题。
针对问题1，我们首先根据附件一提供的景区信息，进一步收集了全国5A景区的经纬度信息，并通过粒子群算法和原点映射绘制了相对分布图和散点图，可视化了这些景区在全国的分布情况，最后得出5A景区的分布与地理位置有关的结论。在此基础上，运用DESTEP分析模型和PESTEL分析模型对影响景区评价水平的因素及其相关性进行了分析。
对于问题2，首先运用层次分析法和熵值法对问题1中可能影响景区评价水平的因素进行量化和赋权，最后得出景区综合得分的公式，构建了衡量景区接待能力的指标体系。最后，采用K-Means聚类分析方法对景区进行分类，最终确定所收集数据中景区的综合得分。
针对问题3，首先利用元胞自动机模型和Mersenne旋转算法，针对不同风险等级的景区提出了一种差分流量限制方案，并进行了景区客流模拟。然后对流量限制模型进行了扩展，在区分流量限制的背景下考虑了疫情爆发时影响游客紧急疏散的各种因素，并通过随机构建二维栅格图对这些因素进行了模拟，比较了不同情况下游客的疏散状况。
针对问题4，首先应用灰狼多目标动态规划算法和NSGA-II多目标优化模型对所构建的定量模型进行分析，然后利用BP卷积神经网络算法对所构建的多目标规划优化模型进行预测检验，最后得出区域内潜在的疫情风险越低，景区的限制越低的结论。区域内潜在的传染病风险越低，流量越小，景区的最大瞬时承载能力越大，从而增加景区的收入，提高游客的旅游体验。
针对问题5，我们首先在前四个问题调研分析的基础上，结合国家文旅办和省、自治区、直辖市政府发布的政策文件传达的思路，提出了以政府为主体，针对不同风险等级景区的差异化管理方案。
对于这一问题的求解，我们建立并应用的模型较好地解决了问题中的五个问题。通过小波分析算法、RBP神经网络算法和GRNN-PNN联合神经网络算法对模型进行误差测试和灵敏度分析，证明所建立的模型具有更可靠的灵敏度和更高的可靠性。关键词：粒子群算法、K均值聚类分析算法、元胞自动机模型、灰狼多目标动态规划算法、NSGA-II多目标优化模型、BP卷积神经网络算法、小波分析算法、RBP神经网络算法、GRNN-PNN联合神经网络算法。

模型假设：

我们假设景区内的游客在游览过程中是在景区工作人员的指挥下，且速度恒定，不会有游客长时间滞留在景区内受阻。
在众多影响景区综合得分的因素中，除文中提到的因素外，其他因素的影响可以忽略不计。

问题分析：

数据分析
应用粒子群算法、DESTEP和PESTEL分析模型、层次分析法、熵法、三维置信椭球法、KMeans聚类算法、元胞自动机算法、Gray Wolf多目标动态规划算法、NSGA-II多目标优化算法、BP卷积神经网络预测算法、小波分析算法、，针对RBP神经网络算法、GRNN神经网络算法、PNN神经网络算法这一问题进行数据分析。算法有：小波分析算法、RBP神经网络算法、GRNN神经网络算法、PNN神经网络算法。

问题一分析
根据附件一提供的景区信息，进一步收集我国5A景区的经纬度信息，并使用粒子。通过聚类算法和原点映射绘制出这些景区在全国的相对分布图和散点图，并利用DESTEP分析模型和PESTEL分析模型分析了影响景区评价水平的因素及其相互关系。

问题二分析
运用AHP层次分析法和熵值法对可能影响景区评价水平的因素进行量化赋权，最终得出景区综合评分公式作为衡量景区接待能力的指标，并进行方差和标准差检验。最后，利用K-Means聚类分析方法对景区进行分类。

问题三分析
在流行病背景下，利用元胞自动机模型和Mason旋转算法，提出了一种针对不同风险等级景区的差分流量限制方案，并进行了仿真。然后将该模型推广到考虑不同流量限制条件下突发疫情下游客紧急疏散的各种影响因素，并通过随机构建二维栅格地图对这些影响因素进行模拟和建模。

问题4的分析
该问题是一个多目标动态规划优化问题。通过引入Gray Wolf多目标规划算法和NSGA-II多目标优化模型对构建的定量模型进行分析，并利用BP卷积神经网络算法对构建的多目标规划优化模型进行预测检验。

问题五分析
在对前四个问题进行调研分析的基础上，结合国家文旅办和省、自治区、直辖市政府政策文件传达的思路，提出了不同风险等级景区的差异化管理方案，以政府为主体。

模型的建立与求解整体论文缩略图

在这里插入图片描述

全部论文及程序请见下方“ 只会建模 QQ名片” 点击QQ名片即可

程序代码：

%% Download Data
data=load('quanguojingqujingweidu.txt');
cityCoor=[data(:,2) data(:,3)]; % city coordinate matrix
figure
plot(cityCoor(:,1),cityCoor(:,2),'ms','LineWidth',2,'MarkerEdgeColor','k','MarkerFace
Color','g')
legend('Scenic location')
ylim([4 78])
title('scenic distribution map','fontsize',12)
xlabel('km','fontsize',12)
ylabel('km','fontsize',12)
%ylim([min(cityCoor(:,2))-1 max(cityCoor(:,2))+1])
grid on
% Calculate inter-city distance
n=size(cityCoor,1); %Number of cities
cityDist=zeros(n,n); %city distance matrix
for i=1:n
for j=1:n
if i~=j
cityDist(i,j)=((cityCoor(i,1)-cityCoor(j,1))^2+...
(cityCoor(i,2)-cityCoor(j,2))^2)^0.5;
end
cityDist(j,i)=cityDist(i,j);
end
end
nMax=200; %evolution number
indiNumber=1000; %Number of individuals
individual=zeros(indiNumber,n);
%^ Initialize particle position
for i=1:indiNumber
individual(i,:)=randperm(n); 
end
%% Calculation of population fitness
indiFit=fitness(individual,cityCoor,cityDist);
[value,index]=min(indiFit);
tourPbest=individual; %Current individual best
tourGbest=individual(index,:) ; %Current global optimum
recordPbest=inf*ones(1,indiNumber); %individual best record
recordGbest=indiFit(index); %Group best record
xnew1=individual;
%% Loop to find the optimal path
L_best=zeros(1,nMax);
for N=1:nMax
N
% Calculate the fitness value
indiFit=fitness(individual,cityCoor,cityDist);
% Update current best and historical best
for i=1:indiNumber
if indiFit(i)<recordPbest(i)
recordPbest(i)=indiFit(i);
tourPbest(i,:)=individual(i,:);
end
if indiFit(i)<recordGbest
recordGbest=indiFit(i);
tourGbest=individual(i,:);
end
end
[value,index]=min(recordPbest);
recordGbest(N)=recordPbest(index);
%% Crossover operation
for i=1:indiNumber
% Crossover with individual optimal
c1=unidrnd(n-1); %generate crossover bits
c2=unidrnd(n-1); %generate crossover bits
while c1==c2
c1=round(rand*(n-2))+1;
c2=round(rand*(n-2))+1;
end
chb1=min(c1,c2);
chb2=max(c1,c2);
cros=tourPbest(i,chb1:chb2);
ncros=size(cros,2); 
% Delete the same element as the intersection area
for j=1:ncros
for k=1:n
if xnew1(i,k)==cros(j)
xnew1(i,k)=0;
for t=1:n-k
temp=xnew1(i,k+t-1);
xnew1(i,k+t-1)=xnew1(i,k+t);
xnew1(i,k+t)=temp;
end
end
end
end
% Insert crossover area
xnew1(i,n-ncros+1:n)=cros;
Accept if the new path length becomes shorter
dist=0;
for j=1:n-1
dist=dist+cityDist(xnew1(i,j),xnew1(i,j+1));
end
dist=dist+cityDist(xnew1(i,1),xnew1(i,n));
if indiFit(i)>dist
individual(i,:)=xnew1(i,:);
end
% Crossover with all optimal
c1=round(rand*(n-2))+1; %generate crossover bits
c2=round(rand*(n-2))+1; %generate crossover bits
while c1==c2
c1=round(rand*(n-2))+1;
c2=round(rand*(n-2))+1;
end
chb1=min(c1,c2);
chb2=max(c1,c2);
cros=tourGbest(chb1:chb2); 
ncros=size(cros,2); 
% Delete the same element as the intersection area
for j=1:ncros
for k=1:n
if xnew1(i,k)==cros(j)
xnew1(i,k)=0;
for t=1:n-k
temp=xnew1(i,k+t-1);
xnew1(i,k+t-1)=xnew1(i,k+t);
xnew1(i,k+t)=temp;
end
end
end
end
% Insert crossover area
xnew1(i,n-ncros+1:n)=cros;
Accept if the new path length becomes shorter
dist=0;
for j=1:n-1
dist=dist+cityDist(xnew1(i,j),xnew1(i,j+1));
end
dist=dist+cityDist(xnew1(i,1),xnew1(i,n));
if indiFit(i)>dist
individual(i,:)=xnew1(i,:);
end
%% Variation operation
c1=round(rand*(n-1))+1; %generate variant bits
c2=round(rand*(n-1))+1; %generate variant bits
while c1==c2
c1=round(rand*(n-2))+1;
c2=round(rand*(n-2))+1;
end
temp=xnew1(i,c1);
xnew1(i,c1)=xnew1(i,c2);
xnew1(i,c2)=temp;
Accept if the new path length becomes shorter
dist=0;
for j=1:n-1
dist=dist+cityDist(xnew1(i,j),xnew1(i,j+1));
end
dist=dist+cityDist(xnew1(i,1),xnew1(i,n));
if indiFit(i)>dist
individual(i,:)=xnew1(i,:);
end
end
[value,index]=min(indiFit);
L_best(N)=indiFit(index);
tourGbest=individual(index,:); 
end
%% Results for graphing
figure
plot(L_best)
title('Algorithm training process')
xlabel('number of iterations')
ylabel('Adaptation value')
grid on
figure
hold on
plot([cityCoor(tourGbest(1),1),cityCoor(tourGbest(n),1)],[cityCoor(tourGbest(1),2),...
cityCoor(tourGbest(n),2)],'ms-
','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','g')
hold on
for i=2:n
plot([cityCoor(tourGbest(i-1),1),cityCoor(tourGbest(i),1)],[cityCoor(tourGbest(i1),2),...
cityCoor(tourGbest(i),2)],'ms-
','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','g')
hold on
end
legend('x')
scatter(cityCoor(:,1),cityCoor(:,2));
title('x','fontsize',10)
xlabel('km','fontsize',10)
ylabel('km','fontsize',10)
grid on
ylim([4 80])

% Transformation of input data into statistical data
close all;clear;
%% Loading data 
path_train = 'TrainData_2021.2.1_2021.2.19.txt';
path_test = 'TestData_2021.2.20_2021.2.27.txt';
[traind_flavor,train_time] = readtxt(path_train);
[testd_flavor,test_time] = readtxt(path_test);
% Prediction of data
% % When the VM is 0 in all moments, the prediction is that these data are also 0
% sum_flavor1 = sum(train_flavor);
% sum_flavor2 = sum(test_flavor);
% ind_zero = find(sum_flavor1==0);
% ind_pre(1:ind_zero) = 0;
% Number prediction for virtual machine specifications that are not 0
data_nuum = length(traind_flavor(:,1));
traind = traind_flavor(1:data_nuum - 7,:);
trainl = traind_flavor(8:data_nuum,:);
% data_nuum2 = length(testd_flavor(:,1));
% testd = testd_flavor(1:data_nuum2 - 6,:);
% testl = testd_flavor(7:data_nuum2,:);
%% Create Network
net = feedforwardnet(31);
net.trainFcn = 'trainbfg';
% net.trainFcn = 'trainlm';
net.trainParam.epochs=1000;%Maximal number of training steps allowed 1000
net.trainParam.max_fail = 10;
% view(net);
%% Training Network
net = train(net,traind',trainl');
%% Test
test_out=sim(net,traind_flavor((data_nuum - 6):data_nuum,:)');
test_out = test_out';
test_out(test_out<10)=0;
test_out = round(test_out);
% % Data Evaluation
N = 15;
for i =1:7
s1(i) = sqrt(sum((sum(test_out(i,:)) - testd_flavor(i,:)). ^2)/N);
s2(i) = sqrt(sum(test_out(i,:). ^2));
s3(i) = sqrt(sum(testd_flavor(i,:). ^2));
score(i) = 1- (s1/(s2+s3));
end
function main()
clc;clear all;close all;
% Use Mexihat function as sample input and output
x=0:0.03:3; %Sample input value
c=2/(sqrt(3). *pi.^(1/4));
d=1/sqrt(2);
u=x/2-1;
targ=d.*c.*exp(-u.^2/2). *(1-u.^2); % sample output value of the objective function
eta=0.02;aerfa=0.735; %Give initial values of network learning rate and momentum 
factor
% Initialize the connection rights wjh the output layer and the hidden layer and the 
connection rights of the hidden layer and the output layer.
% Assume that the number of wavelet function nodes is H, the number of samples is 
P, the number of output nodes is J, and the number of input nodes is I.
H=15;P=2;I=length(x);J=length(targ);
b=rand(H,1);a=rand(H,1); %initialize wavelet parameters
whi=rand(I,H);wjh=rand(H,J); % initialize the weight coefficients.
b1=rand(H,1);b2=rand(J,1); %Threshold initialization;
p=0;
Err_NetOut=[]; % saved errors.
flag=1;count=0;
while flag>0
flag=0;
count=count+1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xhp1=0;
for h=1:H
for i=1:I
xhp1=xhp1+whi(i,h)*x(i);
end
ixhp(h)=xhp1+b1(h);
xhp1=0;
end
for h=1:H
oxhp(h)=fai((ixhp(h)-b(h))/a(h));
end
ixjp1=0;
for j=1:J
for h=1:H
ixjp1=ixjp1+wjh(h,j)*oxhp(h);
end
ixjp(j)=ixjp1+b2(j);
ixjp1=0;
end
for i=1:J
oxjp(i)=fnn(ixjp(i));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wuchayy=1/2*sumsqr(oxjp-targ);
Err_NetOut=[Err_NetOut wuchayy]; %Save the error for each time;
% Solve the wavelet network using BP algorithm, the amount of adjustment of each 
parameter learned each time
for j=1:J
detaj(j)=-(oxjp(j)-targ(j))*oxjp(j)*(1 - oxjp(j));
end
for j=1:J
for h=1:H
detawjh(h,j)=eta*detaj(j)*oxhp(h);
end
end
detab2=eta*detaj;
sum=0;
for h=1:H
for j=1:J
sum=detaj(j)*wjh(h,j)*diffai((ixhp(h)-b(h))/a(h))/a(h)+sum;
end
detah(h)=sum;
sum=0;
end
for h=1:H
for i=1:I
detawhi(i,h)=eta*detah(h)*x(i);
end
end
detab1=eta*detah;
detab=-eta*detah;
for h=1:H
detaa(h)=-eta*detah(h)*((ixhp(h)-b(h))/a(h));
end
% Introduce momentum factor aerfa to speed up convergence and hinder falling into 
local minima.
wjh=wjh+(1+aerfa)*detawjh;
whi=whi+(1+aerfa)*detawhi;
a=a+(1+aerfa)*detaa';
b=b+(1+aerfa)*detab';
b1=b1+(1+aerfa)*detab1';
b2=b2+(1+aerfa)*detab2';
% This algorithm uses sample-by-sample processing instead of data batch processing
p=p+1;
if p~=P
flag=flag+1;
else
if Err_NetOut(end)>0.008
flag=flag+1;
else
figure;
plot(Err_NetOut);
xlabel('number of times the network learned');ylabel('error of network output');
title('network learning error curve','fontsize',20,'color',[0 1 1],'fontname','scribe');
end
end
if count>6000
figure(1);
subplot(1,2,1)
plot(Err_NetOut,'color','b','linestyle','-','linewidth',2.2,...
'marker','^','markersize',3.5);
xlabel('number of times the network learned');ylabel('error of network output');
title('error-curve','fontsize',20,'color',[1 1 1],'fontname','scribe');
subplot(1,2,2)
handle1=plot(x,targ,'color','r','linestyle','--','linewidth',2.2,...
'marker','p','markersize',3.5);
hold on
handle1=plot(x,oxjp,'color','g','linestyle','-.' ,'linewidth',2.2,...
'marker','d','markersize',3.5);
xlabel('sample input value');ylabel('sample target value and network output value');
title('Comparison of target and network output values','fontsize',20,'color',[1 1 
1],'fontname','scribe');
legend('sample target value','network simulation value');
break;
end
end
function y3=diffai(x) % subroutine
y3=-1.75*sin(1.75*x). *exp(-x.^2/2)-cos(1.75*x). *exp(-x.^2/2). *x;
function yl=fai(x) % subroutine
yl=cos(1.75.*x). *exp(-x.^2/2);
function y2=fnn(x) % subroutine
y2=1/(1+exp(-x));