2021年数维杯国际大学生数学建模
A题 新冠肺炎背景下港口资源优化配置策略
原题再现：
2020年初，新型冠状病毒（COVID-19）在全球迅速蔓延。根据世界卫生组织2021年7月31日的报告，新冠病毒疫情对人类的影响可能比原先预期的持续时间更长。在这一流行病的影响下，许多国家遭受了不同程度的经济损失，各行各业也面临着许多前所未有的困难，全球贸易的发展也迎来了许多新的挑战。
作为贸易发展中最重要的环节之一，港口和航运业同样难逃厄运。新冠病毒疫情的爆发暴露出全球港口和航运业在运营、管理方面缺乏系统协调。特别是港口和航运的多维度防疫措施，大大降低了全球港口的运营效率，导致航运服务日益短缺，导致整个行业内发生不可预测的变化。
COVID-19对港口的影响如下：第一，随着亚洲一些国家疫情的恢复，货物订单量持续上升，造成港口严重拥堵。然而，在一些港口极度拥堵的情况下，其他一些港口的货物量变得非常小，出现了港口闲置的现象。这使得港口交通不对称现象更加普遍。第二，新冠病毒还导致世界航运能力不足。例如，波罗的海干散货指数（BDI）在2021年8月20日飙升至4092，这也是过去11年来的最高点。此外，由于港口服务效率低下和航运市场动荡，新冠病毒也加剧了全球供应链的不确定性，导致供应链中断的风险增加。因此，在COVID-19时代，对港口运营效率提出了更高的要求。
在COVID-19疫情持续肆虐的今天，如何降低防控带来的风险成为全球港口航运业可持续管理的关键问题。在这种情况下，我们需要重新设计一个高度灵活的系统和全球网络，以保持利益相关者之间的密切合作关系，包括航运公司、港务局、码头运营商以及物流和港口服务提供商。
为缓解COVID-19影响下港口运营效率低下的问题，请运用数学建模方法解决以下五个问题：
问题1：结合附件1的数据和您的合理假设，请从船公司联合对接的角度出发，提供船公司之间的合作模式，以缓解港口拥堵或过度闲置的问题。你应该深入分析每个码头在每个周期的航运公司。
问题2：请根据附件1的数据，从码头运营商的角度设计船舶最佳时基对接方案，并提出具体的奖惩措施，以有效提高港口的对接效率。
问题3：请从港口泊位、码头运营商和航运公司的角度提供船舶停靠策略。该方案需要在不减少防疫工作时间的前提下，最大限度地缩短船舶在港口的停靠时间。
问题4：能否提供一个充分考虑区域内多个港口合作的策略，从而有效减少港口拥堵。
问题5：如果港口安排1-2天的周末休息，您能否提出一个更通用的港口呼叫优化模型，并使用附件1中的数据进行演示。
整体求解过程概述(摘要)
随着新型冠状病毒（COVID-19）在全球的传播，港口和航运业受到严重打击。主要问题是如何参考附件一给出的船舶载重数据和泊位通过能力，优化港口资源配置，尽可能提高港口运营效率。本文构建了一个异构船舶资源配置优化模型来解决该问题。
对于问题1：首先，从航运公司的角度建立了基于航运公司的所有船舶期望值的表达函数，建立了泊位利用率最大的异质船舶资源配置优化模型，通过考虑所有航运公司的所有船舶的能力和所有泊位在一个周期内的停泊能力，建立了航运公司的最大期望值和最大公平性。并用遗传算法进行求解，结果表明，一个周期内所有泊位利用率之和为5489，各航运公司期望值为3489。
针对问题2：考虑港口航运公司船舶交接时间因素，建立了基于软时间窗的异构船舶资源配置优化模型，设计了基于时间窗的船舶加班惩罚策略。在此基础上，提出了基于插值和聚集的优化策略。结果表明，船舶靠泊计划平均完成率为82%，泊位空置率为14.38%。
针对问题3：在模型1和模型2的基础上，建立了以最小停靠时间为主要目标的基于停靠时间的异构船舶港口泊位资源分配模型。并采用快速非支配排序策略和淘汰策略对模型进行求解。
对于问题4：我们将该问题的维度提升到一个区域内多个港口之间的合作，并考虑多个港口的多个泊位与不同航运公司之间的合作，以改进我们之前构建的模型。采用“三数中立”策略和“淘汰策略”策略对模型进行求解。结果表明，当各端口可用率为65%时，总体效率最大。
对于问题5：基于问题4的模型，我们添加了端口周末休息约束。并采用“三数中立”策略和“淘汰策略”策略对模型进行求解。并进行了仿真实验，验证了基于附录1的模型。最后，我们在“Result presentation”文件夹中展示了结果。
模型假设：
1、本文忽略了泊位类型，所有泊位均为集装箱类型。从文献[1]可知，泊位的水工结构是所述泊位容量的1.5倍，即每个泊位为一艘船舶服务的容量不得超过泊位泊位容量的1.5倍。
2、航运公司为了实现利润最大化，要求每艘船舶只有在满载状态下，即每艘船舶装载到其最大容量时，才能执行航运任务。
3、一船只能选择一个泊位。
4、港口计划运营周期为14天。
5、我们没有计划船的航迹。
6、本文忽略了船舶装卸后泊位的无效等待时间。
7、我方按附件一计算港口泊位工作时间：8:00～18:00。
问题分析：
数据分析
附件1提供了1127艘船舶的代码、载客量、装运周期、平均停靠时间、最近十次停靠时间的差异、最近十次到达时间以及217个泊位的停泊能力。通过对运力数据的整理和分析，发现所有船舶的航运周期分布在1~14天之间，即周期最短的船舶每天到达港口，周期最长的船舶需要14天到达港口，我们发现港口每天早上8:00到晚上18:00接收船舶靠泊，即船舶只能在规定时间内到达港口。
问题1分析
解决这一问题需要从航运公司的角度考虑，既要满足所有航运公司的需要，又要考虑泊位的利用率。船舶的需求因航运公司而异，解决问题的关键是如何综合这些因素进行优化。
问题2的分析
这一问题的解决需要从码头运营商的角度进行充分考虑，在满足泊位利用率的前提下，对船公司的每艘船舶进行时间和泊位规划，因此不同船公司的不同船舶容量和泊位通过能力是需要重点关注的关键因素。解决问题2的关键是将问题1的模型视角转化为终端操作员的模型视角，优化时间进度的处理以满足问题2的要求，最终解决问题。
问题3的分析
该问题侧重于需要同时从码头运营商和航运公司的角度出发，考虑到防疫工作的前提，同时满足不同航运公司的需求，尽可能减少船舶总停靠时间，提出的策略需要将问题1和问题2一并考虑，并添加时间窗等约束条件，通过优化算法求解。
问题4的分析
在问题4中，问题的维度上升到一个区域内多个港口之间的合作，不再局限于单个港口的多个泊位，而是需要考虑多个港口的多个泊位与不同航运公司之间的合作调度，基于博弈论的思想，以多港口合作为关键，以船舶停靠时间为目标，优化船舶停靠时间。
问题5的分析
对问题5的审议将继续深入实际情况，安排港口每个周末休息1-2天。然后在问题4的基础上对模型进行改进，在不延误航运公司船舶需求的前提下，规划各港口的休息时间以满足基本要求，如何解决这一问题是问题的关键。
模型的建立与求解整体论文缩略图
全部论文及程序请见下方“ 只会建模 QQ名片” 点击QQ名片即可
程序代码：
clear all;
close all;
clc;
A=[];
B=[];
C=[];
for gg=1:50
run question1
z=30;
c=184;
Xs1=rand*(3503-3450)+3450;
Xx1=50;
NP=z;
L=c;
Pc=0.8;
Pm=0.1;
G=100;
f=randi([0,1],NP,L)
for k=1:G
%Decode binary to decimal within the domain
for i=1:NP
U=f(i,:);
m=0;
for j=1:L
m=U(j)*2^(j-1)+m;
end
x(i)=Xx1+m*(Xs1-Xx1)/(2^L-1);
Fit(i)= fitness1(x(i));
end
maxFit=max(Fit);
minFit=min(Fit);
rr=find(Fit==maxFit);
fBest=f(rr(1,1),:);
xBest=x(rr(1,1));
Fit=(Fit-minFit)/(maxFit-minFit);
sum_Fit=sum(Fit);
fitvalue=Fit./sum_Fit;
fitvalue=cumsum(fitvalue);
ms=sort(rand(NP,1));
fiti=1;
newi=1;
while newi<=NP
if (ms(newi))<fitvalue(fiti)
nf(newi,:)=f(fiti,:);
newi=newi+1;
else
fiti=fiti+1;
end
end
for i=1:2:NP
p=rand;
if p<Pc
q=randi([0,1],1,L);
for j=1:L
if q(j)==1;
temp=nf(i+1,j);
nf(i+1,j)=nf(i,j);
nf(i,j)=temp;
end
end
end
end
i=1;
while i<=round(NP*Pm)
h=randi([1,NP],1,1);
for j=1:round(L*Pm)
g=randi([1,L],1,1);
nf(h,g)=~nf(h,g);
end
i=i+1;
end
f=nf;
f(1,:)=fBest;
trace1(k)=maxFit;
end
%Second optimization
xBest;
z=30;
c=184;
Xs1=1.50*(rand*(3503
-3450)+3450);
Xx1=50*0.5;
NP=z;
L=c;
Pc=0.8;
Pm=0.1;
G=100;
f=randi([0,1],NP,L)
for kk=1:G
for i=1:NP
U=f(i,:);
m=0;
for j=1:L
m=U(j)*2^(j
-1)+m;
end
x(i)=Xx1+m*(Xs1
-Xx1)/(2^L
-1);
Fit(i)= fitness1(x(i));
end
maxFit=max(Fit);
minFit=min(Fit);
rr=find(Fit==maxFit);
fBest=f(rr(1,1),:);
xBest=x(rr(1,1));
Fit=(Fit
-minFit)/(maxFit
-minFit);
sum_Fit=sum(Fit);
fitvalue=Fit./sum_Fit;
fitvalue=cumsum(fitvalue);
ms=sort(rand(NP,1));
fiti=1;
newi=1;
while newi<=NP
if (ms(newi))<fitvalue(fiti)
nf(newi,:)=f(fiti,:);
newi=newi+1;
else
fiti=fiti+1;
end
end
for i=1:2:NP
p=rand;
if p<Pc
q=randi([0,1],1,L);
for j=1:L
if q(j)==1;
temp=nf(i+1,j);
nf(i+1,j)=nf(i,j);
nf(i,j)=temp;
end
end
end
end
i=1;
while i<=round(NP*Pm)
h=randi([1,NP],1,1);
for j=1:round(L*Pm)
g=randi([1,L],1,1);
nf(h,g)=~nf(h,g);
end
i=i+1;
end
f=nf;
f(1,:)=fBest;
trace2(kk)=maxFit;
end
xBest;
for h=1:100
trace(h)=trace1(h)*trace2(h)*variance_sum;
end
figure
plot(trace)
xlabel('Number of iterations')
ylabel('Objective function value')
title('Evolution curve of fitness value of objective function')
A=[A,trace(100)];
% saveas(gcf,['Code\Question1\','fig',num2str(gg)])
end
B=sort(A);
figure
plot(B)
[~,a]=find(A == min(A))
[~,b]=find(A == max(A))
c=sum(A)/50;
plot(B)