一、模型建立

符号定义：
g——港口位置 $g$ $\in$ $G,1,2,\dots n$
$i$——泊位位置$i$ $\in$ $I,1,2,\dots n$
$j$——船舶位置$j$ $\in$ $J,1,2,\dots n$
$Twj$——船j的等待时间
$Tkj$——船j的卸货时间（一个集装箱的卸货时间是15分钟）
$Taj$——船j的到达时间
$Tlj$——船j的离开时间（ $Tlj=Taj+Tkj+Twj$）
$Dj$——船j吃水深度
$Lj$——船j长度
$Wgi$——港口g泊位i的水深
$Lgi$——港口g中泊位i的长度
$X{{}_g^p}_{ij}$=1船j在港口g的泊位i中第p接受服务， =0船j在港口g的泊位i中不是第p接受服务
$C_{sj}$——船j的航行成本
$C_{wj}$——船j的等待成本（8000元/小时）
$C_{tj}$——船j到达港口后的集装箱运输成本
$N_j$——船j的集装箱总量（z:重箱量； k:空箱量）
$P$——船舶单位航行成本（5元/（t*km））
$A_i$——目的港口的坐标 $A_i^{{}^{\prime }}$——空箱到达目的港口的坐标
$C_{SZ}$——重箱运输成本（单位成本15元/km）
$C_{sk}$——空箱运输成本（单位成本10元/km）
$d_j$——船j重箱商品实际要达到的地理坐标
P——船舶靠港次序（Z——船舶靠港总次序）
船舶经济航速为40km/h

数学模型如下：

（1）$\min T=\sum _{g\in G}\sum _{i\in I}\sum _{j\in J}\sum _{p\in Z}\left(T_{wj}+T_{kj}\right)X_{g_{ij}}^p$ 使船舶在港时间最短
（2）$\min C=\sum _j\left(C_{sj}+C_{wj}+C_{tj}\right)$ 船公司的货运成本最小
（3）$T_{wj}=T_{lj}-T_{aj}-T_{kj}\ge 0$ 船舶等待时间大于等于0
（4）$X_{g_{ij}}^p\times X_{g_{ij}}^{p+1}\times (T_{lj}-T_{aj})\ge 0$ 必须要求前一艘离开另一艘才开始服务
（5）$C_{sj}=N_j\times p\times A_j$ 船 j的航行成本=集装箱总数单位航行距离成本距离
（6）$C_{wj}=C_t\times T_{wj}$ 船j的等待成本=单位等待成本*等待时间
（7）$N=K+Z$ 集装箱总量=重箱数+空箱数
（8）停泊船的吃水深度小于泊位水深；船的长度小于泊位长度，即船都有靠泊的基本条件。
（1）（2）为目标函数，（3）——（8）为约束条件

其实这个问题非常简单，就是一个指派问题，给每一艘船指定一个港口，然后就可以计算得到每个船的等待时间、航行距离等等信息。知道了这个问题，程序就非常好解决啦。

二、模型数据

营口港与大连港的距离为210千米，为简化作出以下示意图

2个港口的情况（在这一部分中，为了计算方便将原数据按一定比例缩小）

船舶数据

三、模拟退火算法编写

主文件：main.m

clc;
clear;
close all;
%% 算法参数
T0=5000;   % 初始温度
Tend=1e-3;  % 终止温度
L=500;    % 各温度下的迭代次数（链长）
q=0.9;    %降温速率

num_bowei = [3 3];%泊位数量
num_boat = 10; %船的数量
N = num_boat;

%% 初始解
S1 = ceil(sum(num_bowei)*rand(1,num_boat));

%% 输出随机解的路径和总距离

Rlength=PathLength(S1);
disp(['总距离：',num2str(Rlength)]);

%% 计算迭代的次数Time
% Time=ceil(double(solve(['1000*(0.9)^x=',num2str(Tend)])));
Time=132;
count=0;        %迭代计数
Obj=zeros(Time,1);         %目标值矩阵初始化
track=zeros(Time,N);       %每代的最优路线矩阵初始化
%% 迭代
while T0>Tend
    count=count+1;     %更新迭代次数
    temp=zeros(L,N+1);
    for k=1:L
        %% 产生新解
        S2=NewAnswer(S1);
        %% Metropolis法则判断是否接受新解
        [S1,R]=Metropolis(S1,S2,T0);  %Metropolis 抽样算法
        temp(k,:)=[S1 R];          %记录下一路线的及其路程
    end
    %% 记录每次迭代过程的最优路线
    [d0,index]=min(temp(:,end)); %找出当前温度下最优路线
    if count==1 || d0<Obj(count-1)
        Obj(count)=d0;           %如果当前温度下最优路程小于上一路程则记录当前路程
    else
        Obj(count)=Obj(count-1);%如果当前温度下最优路程大于上一路程则记录上一路程
    end
    track(count,:)=temp(index,1:end-1);  %记录当前温度的最优路线
    T0=q*T0;     %降温
%     fprintf(1,'%d\n',count)  %输出当前迭代次数
end
%% 优化过程迭代图
figure
plot(1:count,Obj)
xlabel('迭代次数')
ylabel('距离')
title('优化过程')

Result(track(end,:));

目标函数计算

function len=PathLength(Chrom)
%% 计算各个体的路径长度
% 输入：

load = [
    150
    100
    200
    100
    150
    150
    200
    200
    200
    100];

num_zhong = [120
50
100
70
100
110
150
160
160
70
];
num_kong = [30
50
100
30
50
40
50
140
40
30
];


arrive_time = 6:1:15;
num_bowei = [3 3];%泊位数量
num_boat = 10; %穿的数量

%仿真求解适应度值
[NIND,num_boat]=size(Chrom);

wait_time = zeros(1,num_boat);
departure_time = zeros(1,num_boat);
service_time = zeros(1,sum(num_bowei));

c_wait = zeros(1,num_boat); %等待费用
c_sail = zeros(1,num_boat); %航行成本
c_transfer = zeros(1,num_boat); %转移成本


x=Chrom;
for j=1:1:num_boat
    if x(j) > num_bowei(1) % 选择第二个港口，到达时间要延迟
        arrive_time(j) = arrive_time(j) + 210/40;
        c_transfer(j) = (num_zhong(j) * 15 + num_kong(j) * 10) * 210;
        c_sail(j) = load(j) * 5 * 210;
    end
    if arrive_time(j) > service_time(x(j)) % 没有被占用,直接服务
        wait_time(j) = 0;
        departure_time(j) = arrive_time(j) + load(j)*15/60;
        service_time(x(j)) = departure_time(j);
    else % 已经被占用了
        wait_time(j) = service_time(x(j)) - arrive_time(j);
        departure_time(j) = arrive_time(j) + wait_time(j) + load(j)*15/60;
        service_time(x(j)) = departure_time(j);
    end
    
end
c_wait = wait_time .* 8000;
%len = sum(wait_time);
len = sum(c_wait) + sum(c_sail) + sum(c_transfer);

四、运行结果

非常简洁。

最后，博主专注于论文的复现工作，有兴趣的同学可以私信共同探讨。相关代码已经上传到资源共享，点击我的空间查看分享代码。