動車組檢修問題 數學建模

前言

暑期數學建模已經接近尾聲，淘汰賽進行到第二輪，這次的題目是《動車組檢修問題》。題目一發出來后，按照慣例各種找資源、找文獻。最后鎖定了題目來源 同濟大學校內數模競賽2019年 B題 , 仔細對比發現，題目在到站時間上略有區別。

題目展示

動車組運用所是負責對動車進行檢修、養護等工作的場所，全國已建成的動車 運用所已超過 50 個。動車組的檢修根據行駛情況被划分成不同檢修等級，不同 等級對應不同的工序。

問題一

如圖 1 所示，動車組的一次檢修包括 a，b，c 三個工序。每個工序擁有的作業車 間和需要花費的時間如表 1 所示，相同工序不同車間的耗費時間相同。動車組按 a→b→c 順序進行檢修，完成一個檢修工序后駛入下一個有空閑位置的車間進行 下一個檢修工序，若下一個工序所有車間都處於占用狀態，則動車組需要在上一 個車間中等待。動車運用所某 12 小時內每十五分鍾來 1 輛待檢修的動車，按照 目前的車間設置，維修完所有這些動車組總共需要多長時間？請給出安排檢修的 最佳方案。假設第一輛動車組抵達動車運用所時，所有檢修車間都是空閑的，且 車間之間的轉換時間忽略不計。

表1 檢修基礎數據

工序類別	a	b	c
車間數量（個）	3	8	5
耗費時間（小時）	1	2	1.5

問題二

事實上，如表2所示，不同類型的動車組每個工序需要花費的時間是不一樣的。請根據附件中附表一到達動車運用所的動車信息，計算維修完這些動車的總時間。

表2 檢修基礎數據

動車類別 \ 工序類別	a	b	c
CRH2	1	2	1.5
CRH3	0.8	2.4	0.5
CRH5	1.3	2.5	1.5
CRH6	1	2.7	0.3

附表一

到站時間	動車類別
00:16	CRH2
00:47	CRH5
01:22	CRH2
02:00	CRH6
02:21	CRH3
03:02	CRH6
03:31	CRH2
03:59	CRH5
04:04	CRH3
04:27	CRH3
05:09	CRH6

問題三

根據列車的行駛時間、歷程和檢修周期，動車組的檢修被划分成不同檢修等級 I~V，如表 3 不同的檢修等級對應不同的工序組合。工序 d 與 e 分別設有車間 3 和 2 個，相同工序不同車間的耗費時間相同。表 4 為不同動車類別的每個工序需 要耗費的時間。根據附表二的到所列車信息，計算檢修完這些列車需要的總時間 是多少？

表3 檢修等級

檢修等級	對應工序組合
Ⅰ	a → b
Ⅱ	a → b → c
Ⅲ	a → b → d
Ⅳ	a → c → d → e
Ⅴ	a → b → c → d → e

表4 檢修基礎數據

動車類別 \ 工序類別	a	b	c	d	e
CRH2	1	2	1.5	4	7
CRH3	0.8	2.4	0.5	4.8	6.5
CRH5	1.3	2.5	1.5	3	6
CRH6	1	2.7	0.3	5	7

附表二

到站時間	動車類別	檢修類型
0:16	CRH2	Ⅳ
0:47	CRH5	Ⅱ
1:22	CRH2	Ⅱ
2:00	CRH6	Ⅰ
2:20	CRH3	Ⅲ
3:05	CRH6	Ⅱ
3:31	CRH2	Ⅴ

解題思路

解題目標

根據題目所述場景，列車按照到站時刻表依次進入到 A 工序各車間進行檢修。若 A 工序車間滿了，則需要等待。檢修完后，列車又會先后離開 A 工序的各個車間。然后再先后抵達 B 工序各個車間進行檢修，然后先后離開B 車間，再先后進入 C 工序車間，最后先后離開C 工序車間。而C工序車間里最后一輛車離開車間時，標志着檢修任務完成。我們就是要讓最后一輛車的離開時間盡可能地早，縮短檢修總耗時。

從進出站時間入手

仔細觀察並思考，就會發現，其實對於每個工序來說，都存在列車進站時間和出站時間。對於A 工序來說，進站時間由題目所給到站時間表決定，這是已知條件。而A 工序的出站時間，是由檢修調度安排決定的，但無論怎樣，列車都會檢修完立馬離開工序車間。於是這里有會有一個出站時間。對於B 工序車間來說，先離開A 工序車間的列車，必然先到達B工序車間。A工序車間的出站時間就可以等同於B工序車間的進站時間（這里可以先不管車間，只考慮工序，因為工序是總的進出）。同理，B 的出站時間也是C的進站時間。最后C的出站時間也將計算出來，這樣一環套一環地去計算，題目就解出來了。

確定編程框架

不難發現，通過上面的分析，我們將三個工序都抽象成了同一個對象，那就是，給定一個進站時間，安排檢修排隊策略，最后輸出出站時間。這就是總體的編程的思路！

各工序車間檢修安排策略

根據上面的流程圖，每個車間都需要依據進站時間安排檢修，當車間有空閑時，可以安排檢修；而當車間無空閑時，則需要等待，直至有空閑的車間出現，如此循環往復，直到所有列車全部檢修完畢，該工序就可以輸出出站時間了。這里有這么幾個問題需要用代碼實現判斷:

• 如何確定當前工序的車間是否有空閑？
• 如何表示列車等待過程，時間怎么計算？
• 如何計算列車的出站時間？
• 列車在哪個車間維修，怎么記錄？

以上問題，都是檢修安排策略的主要問題。我們來整理下編程思路：

首先，列車是按時序來進站的，如何來表示時序？沒錯，給定一個數組，數組元素編號就可以表示是第幾分鍾了（這里的時間單位可以是分鍾也可以是刻鍾，即15分鍾，問題一中就是15分鍾為一個單位考慮，問題二三則以1分鍾為單位考慮）。同理，車間也可以用這樣一個數組表示。考慮問題時，這二者是要結合在一起的。因為在某個時刻，1車間占用，那么就需要考慮去2車間，若也占用，那就去3車間，再找不到，那就只好去下一個時刻了，依舊判斷1車間還在占用沒，若占用，去2車間，... 總會找到某個時刻某個車間是空閑的，此時就安排檢修。於是該車間的后續幾個時刻都將會被占用。這里可以用一個二維數組來記錄上述情形。

給一個M*3的二維矩陣，這里M足夠大，因為由於等待，M不知道要要多大，所以盡量大一些，比如1000（15分鍾為單位）。

第一輛車來的時候，如下圖所示：

第一個15分鍾，車間全是0，表示全是空閑，於是就安排上了，放在第一個車間檢修。

第二輛車來的時候：

此時，第二個15分鍾時，第一趟車已經在檢修，所以占了連續15*4=60分鍾的車間使用時長，於是只好考慮去二車間。

同理，畫出第三輛車的情況：

到了第四輛車來的時候，就出現問題了，此時三個車間全滿，只好去下一個時刻找空位（圖中時間占用沒畫完整）：

這就是排隊等待策略！以上過程通過循環來寫。判斷每一個時刻下，該行為0的所有列，然后將找到的第一個0，作為安排檢修的對象，將該列后續的維修時長個數的值全置為1，表示占用。當找不到為0的列時，則自動跳到下一個時刻繼續找，直到找到為止。然后下一輛車繼續，...... 需要注意的是，起始時刻是由列車的進站時間來決定的，也就是說從第幾個行開始查找。這里可以考慮將列車的進站時間寫成一個數組，按編號讀取到的時間也就是第幾輛車的進站時間。此外，這里找空閑車間的策略是以最小號為選擇對象，即找到的不為1的序列中最小的列號對應車間為當前安排的車間序號。

一直循環下去當然得有一個停止條件，那就是車輛全部檢修完畢，車間空閑時。因此還需要寫一個記錄檢修完成計數矩陣。

總的來說，編程的思路是這樣的流程：

最后一定要將出站時間升序排列！！因為原則是先到先檢修，時間最短，有可能你先進站，結果你后出站，下一工序時，你就不是最前面檢修的對象了。做升序排列目的是讓時間重新恢復遞增狀態，這樣才能完成兩個工序之間的交接了。

以上算法采取的是每一步都按照最優的策略去安排檢修，是一種貪婪算法。整體算法流程如下：

問題解答

有了以上分析后，其實問題一二三都是同一類問題。

• 對於問題一，列車進站時間均勻間隔，每15分鍾一趟，每輛車在同一工序的耗時一樣，屬於最特殊的一種情形。
• 對於問題二，列車進站時間不再均勻，且耗時不一樣，但這並不影響模型的整體框架，我們只需要將車次和維修時間綁定在一起，來什么車就按照什么車的時間安排占位長度。
• 對於問題三，列車進站時間不同，耗時不一樣，車間增加，維修等級還決定了什么車間去，什么車間不去。看似十分復雜，但是稍加修改就可以繼續使用所建立的模型。將車次、耗時、維修等級全部綁定在一起，如果按照等級不去某個車間，則將車間耗時記作0即可，照樣可以得到各車間，各個時刻，各個車次的維修耗時。

問題一答案

甘特圖

列車進出站時間表太長了，在這里就不貼了。

檢修耗時

根據模型計算，考慮從0點列車進站，到最后一輛車出站，全部檢修完需要 (80*15-0)/60 = 20 個小時，屆時已經是20:00。

問題二答案

甘特圖

列車進出站時間表

列車編號	A工序進站時間	A工序出站時間	B工序進站時間	B工序出站時間	C工序進站時間	C工序出站時間
R1	0:16	1:16	1:16	3:16	3:16	4:46
R2	0:47	2:05	2:05	4:35	4:35	6:05
R3	1:22	2:22	2:22	4:22	4:22	5:52
R4	2:00	3:00	3:00	5:42	5:42	6:00
R5	2:21	3:09	3:09	5:33	5:33	6:03
R6	3:02	4:02	4:02	6:44	6:44	7:02
R7	3:31	4:31	4:31	6:31	6:31	8:01
R8	3:59	5:17	5:17	7:47	7:47	9:17
R9	4:04	4:52	4:52	7:16	7:16	7:46
R10	4:32	5:20	5:20	7:44	7:44	8:14
R11	5:09	6:09	6:09	8:51	8:51	9:09

檢修耗時

根據模型計算，考慮從00:16列車進站，到最后一輛車出站，全部檢修完需要 557 - 16 = 541 分鍾，屆時已經是09:17。

問題三答案

甘特圖

列車進出站時刻表

列車編號	A工序進站時間	A工序出站時間	B工序進站時間	B工序出站時間	C工序進站時間	C工序出站時間	D工序進站時間	D工序出站時間	E工序進站時間	E工序出站時間
R1	0:16	1:16	1:16	1:16	1:16	2:46	2:46	6:46	6:46	13:46
R2	0:47	2:05	2:05	4:35	4:35	6:05	6:05	6:05	6:05	6:05
R3	1:22	2:22	2:22	4:22	4:22	5:52	5:52	5:52	5:52	5:52
R4	2:00	3:00	3:00	5:42	5:42	5:42	5:42	5:42	5:42	5:42
R5	2:20	3:08	3:08	5:32	5:32	5:32	5:32	10:20	10:20	10:20
R6	3:05	4:05	4:05	6:47	6:47	7:05	7:05	7:05	7:05	7:05
R7	3:31	4:31	4:31	6:31	6:31	8:01	8:01	12:01	12:01	19:01

檢修耗時

根據模型計算，考慮從00:16列車進站，到最后一輛車出站，全部檢修完需要 1141 - 16 = 1125分鍾，屆時已經是19:01。

源碼展示

貪婪算法模型代碼

function df_A  = mygreedy3(time,n_workshop_A,time_cost_A,n_car,in_ind_A)
time_cost_A = time_cost_A(in_ind_A);
M = 10000;
is_finish_A = zeros(1,n_car); % 記錄所有動車的工藝完成狀態，1為完成，0為未完成
is_work_A = zeros(M,n_workshop_A); % 記錄各個時刻車間的占用情況
input_time_record_A = []; % 記錄動車進站時間
output_time_record_A = []; % 記錄動車出戰時間
workshop_ind_record_A = []; % 記錄檢修車間編號
n_arrival_car_A = 0; % 記錄檢修完的車輛數

for time_A = 1:length(time)
    num_work_A = sum(is_work_A(time(time_A),:)); %計算time時刻的車間運作數量
    num_finish_A = sum(is_finish_A); % 計算time時刻已經發車的數量
    %車間全部空閑 且 所有車已經完成檢修
    if  n_arrival_car_A == n_car
        break
        %車間空閑 且 所有車未完成檢修
    elseif num_work_A < n_workshop_A && num_finish_A < n_car
        [~,col] = find(is_work_A(time(time_A),:) == 0);
        free_id = col(1);
        is_work_A(time(time_A):time(time_A)+time_cost_A(time_A),free_id) = 1;
        is_finish_A(1,n_arrival_car_A+1) = 1;
        n_arrival_car_A = n_arrival_car_A + 1;
        output_time_record_A = [output_time_record_A,time(time_A) + time_cost_A(time_A)];
        input_time_record_A = [input_time_record_A,time(time_A)];
        workshop_ind_record_A = [ workshop_ind_record_A ,free_id];
        %車間無空閑 且 所有車未完成檢修
    elseif num_work_A == n_workshop_A && num_finish_A < n_car
        [~,col] = find(is_work_A(time(time_A),:) == 0);
        k=0;
        while isempty(col)
            k=k+1;
            [~,col] = find(is_work_A(time(time_A)+k,:) == 0);
        end
        free_id = col(1);
        is_work_A(time(time_A)+k:time(time_A)+k+time_cost_A(time_A),free_id) = 1;
        is_finish_A(1,n_arrival_car_A+1) = 1;
        n_arrival_car_A = n_arrival_car_A + 1;
        output_time_record_A = [output_time_record_A,time(time_A)+k + time_cost_A(time_A)];
        input_time_record_A = [input_time_record_A,time(time_A)+k];
        workshop_ind_record_A = [ workshop_ind_record_A ,free_id];
    end
end
% 將結果存到結構體中
for i = 1:length(input_time_record_A)
    df_A(i).input_time = input_time_record_A(i);
    df_A(i).output_time = output_time_record_A(i);
    df_A(i).repair_time =output_time_record_A(i) - input_time_record_A(i);
    df_A(i).workshop_ind = workshop_ind_record_A(i);
    df_A(i).out_ind  = in_ind_A(i);
end
%數據按出站時間升序排列
[~,sort_ind] = sort([df_A.output_time]);
df_A = df_A(sort_ind);
end

問題三代碼

由於問題 一、二、三具有相似性，而問題三最具復雜性，所以給出問題三的代碼，其他兩問稍作精簡和修改即可使用，代碼如下：

%% 准備存儲空間
clc , clear , close all;

%% 變量初始化

[~,~,time ]= xlsread('列車信息.xlsx','列車時刻表');
[~,~,time_cost ]= xlsread('列車信息.xlsx','列車檢修耗時表');
[~,~,rank ]= xlsread('列車信息.xlsx','檢修等級');

train_time = cell2mat(time(2:end,2)); % 列車發車時間
train_id = cell2mat(time(2:end,4)); % 列車編號
train_rank = cell2mat(time(2:end,6));% 列車檢修類型
time_cost = cell2mat(time_cost(2:end,2:end))*60; % 列車檢修耗時
rank = cell2mat(rank(2:end,2:end)); % 檢修類型檢修路徑


% 各工序檢修耗時（單位：1分鍾）

time_cost_A = time_cost(train_id,1).*rank(train_rank,1);
time_cost_B = time_cost(train_id,2).*rank(train_rank,2);
time_cost_C = time_cost(train_id,3).*rank(train_rank,3);
time_cost_D = time_cost(train_id,4).*rank(train_rank,4);
time_cost_E = time_cost(train_id,5).*rank(train_rank,5);


n_car = length(train_time); % 計算檢修車輛數

% 各工序車間數
n_workshop_A = 3;
n_workshop_B = 8;
n_workshop_C = 5;
n_workshop_D = 3;
n_workshop_E = 2;

%% 計算 A 車間的檢修情況
df_A = mygreedy3(train_time,n_workshop_A,time_cost_A,n_car,1:n_car);

%% 計算 B 車間的檢修情況
df_B = mygreedy3(cat(1,df_A.output_time),n_workshop_B,time_cost_B,n_car,cat(1,df_A.out_ind));

%% 計算 C 車間的檢修情況
df_C = mygreedy3(cat(1,df_B.output_time),n_workshop_C,time_cost_C,n_car,cat(1,df_B.out_ind));

%% 計算 D 車間的檢修情況
df_D = mygreedy3(cat(1,df_C.output_time),n_workshop_D,time_cost_D,n_car,cat(1,df_C.out_ind));

%% 計算 E 車間的檢修情況
df_E = mygreedy3(cat(1,df_D.output_time),n_workshop_E,time_cost_E,n_car,cat(1,df_D.out_ind));

%% 整理各列車進出時間
clear R
for i = 1:n_car
    R(i).A= df_A(find([df_A.out_ind]==i));
    R(i).B= df_B(find([df_B.out_ind]==i));
    R(i).C= df_C(find([df_C.out_ind]==i));
    R(i).D= df_D(find([df_D.out_ind]==i));
    R(i).E= df_E(find([df_E.out_ind]==i));
end

%% 結果輸出
time_max = max (cat(1,df_E.output_time));
time_min = min (cat(1,df_A.input_time));
fprintf('列車全部檢修完畢需要%d分鍾，屆時為  %s\n',time_max - time_min,min2time(time_max));

%% 計算列車的進出站時刻表
for i = 1:n_car+1
    if i == 1
        train_io{i,1} = '列車編號';
        train_io{i,2} = 'A工序進站時間';
        train_io{i,3} = 'A工序出站時間';
        train_io{i,4} = 'B工序進站時間';
        train_io{i,5} = 'B工序出站時間';
        train_io{i,6} = 'C工序進站時間';
        train_io{i,7} = 'C工序出站時間';
        train_io{i,8} = 'D工序進站時間';
        train_io{i,9} = 'D工序出站時間';
        train_io{i,10} = 'E工序進站時間';
        train_io{i,11} = 'E工序出站時間';
    else
        train_io{i,1}= sprintf('R%d',i-1);
        train_io{i,2} =min2time(R(i-1).A.input_time);
        train_io{i,3} =min2time(R(i-1).A.output_time);
        train_io{i,4} =min2time(R(i-1).B.input_time);
        train_io{i,5} =min2time(R(i-1).B.output_time);
        train_io{i,6} =min2time(R(i-1).C.input_time);
        train_io{i,7} =min2time(R(i-1).C.output_time);
        train_io{i,8} =min2time(R(i-1).D.input_time);
        train_io{i,9} =min2time(R(i-1).D.output_time);
        train_io{i,10} =min2time(R(i-1).E.input_time);
        train_io{i,11} =min2time(R(i-1).E.output_time);
    end
end
xlswrite('列車進出站時刻表',train_io,'Sheet1')
%% 畫甘特圖
figure;
set(gcf,'outerposition',get(0,'screensize'));
for i = 1:n_car
    rectangle('position',[(R(i).A.input_time) , (R(i).A.workshop_ind-1), (R(i).A.repair_time) ,1]...
        ,'linewidth',0.5,'facecolor',[146,208,80]/255)
    text( (R(i).A.input_time+R(i).A.repair_time/4) , (R(i).A.workshop_ind-1+0.5) , sprintf('R%d-A%d',i,R(i).A.workshop_ind))
    
    rectangle('position',[(R(i).B.input_time) , (R(i).B.workshop_ind + 3 -1), ...
        (R(i).B.repair_time) , 1 ],'linewidth',0.5,'facecolor',[0,176,240]/255)
    text( (R(i).B.input_time+R(i).B.repair_time/4) , (R(i).B.workshop_ind-1 + 3 +0.5) , sprintf('R%d-B%d',i,R(i).B.workshop_ind))
    
    rectangle('position',[(R(i).C.input_time) , (R(i).C.workshop_ind + 3+ 5  -1) , ...
        (R(i).C.repair_time) , 1 ],'linewidth',0.5,'facecolor',[255,192,0]/255)
    text((R(i).C.input_time+R(i).C.repair_time/4) , (R(i).C.workshop_ind-1 + 3 + 5  +0.5) ,sprintf('R%d-C%d',i,R(i).C.workshop_ind))
    
    rectangle('position',[(R(i).D.input_time) , (R(i).D.workshop_ind + 3+5+3  -1) , ...
        (R(i).D.repair_time) , 1 ],'linewidth',0.5,'facecolor',[198,89,17]/255)
    text((R(i).D.input_time+R(i).D.repair_time/4) , (R(i).D.workshop_ind-1 + +3+5+3  +0.5) ,sprintf('R%d-D%d',i,R(i).D.workshop_ind))
    
    rectangle('position',[(R(i).E.input_time) , (R(i).E.workshop_ind + 3+5+3+3 -1 ) , ...
        (R(i).E.repair_time) , 1 ],'linewidth',0.5,'facecolor',[112,48,160]/255)
    text((R(i).E.input_time+R(i).E.repair_time/4) , (R(i).E.workshop_ind-1 + 3+5+3+3  +0.5) ,sprintf('R%d-E%d',i,R(i).E.workshop_ind))
    
end
workname = {'',...
    '工序A-車間1','工序A-車間2','工序A-車間3',...
    '工序B-車間1','工序B-車間2','工序B-車間3','工序B-車間4','工序B-車間5',...
    '工序C-車間1','工序C-車間2','工序C-車間3',...
    '工序D-車間1','工序D-車間2','工序D-車間3',...
    '工序E-車間1','工序E-車間2','工序E-車間3'
    };
xlabel('時間（單位：1分鍾）')
set(gca,'xtick',0:60:1200,'ytick',[0,0.5:16.5],'yticklabel',workname,'fontsize',24)

saveas(gcf,'問題三列車檢修甘特圖.tif')

幾個表格的數據是這樣的：

列車時刻表

到站時間	數值時間	動車類別	列車編號	檢修類型	檢修等級
0:16	16	CRH2	1	IV	4
0:47	47	CRH5	3	II	2
1:22	82	CRH2	1	II	2
2:00	120	CRH6	4	I	1
2:20	140	CRH3	2	III	3
3:05	185	CRH6	4	II	2
3:31	211	CRH2	1	V	5

列車檢修耗時表

	a	c	b	d	e
CRH2	1	1.5	2	4	7
CRH3	0.8	0.5	2.4	4.8	6.5
CRH5	1.3	1.5	2.5	3	6
CRH6	1	0.3	2.7	5	7

檢修等級

檢修等級	工序A	工序B	工序C	工序D	工序E
I	1	1	0	0	0
II	1	1	1	0	0
III	1	1	0	1	0
IV	1	0	1	1	1
V	1	1	1	1	1

分鍾轉時間函數

function time = min2time(minute)
hour = fix(minute/60);
minute = mod(minute,60);
if hour < 10 && minute <10
    time = sprintf('0%d:0%d',hour,minute);
elseif hour >= 10 && minute <10
    time = sprintf('%d:0%d',hour,minute);
elseif hour >= 10 && minute >=10
    time = sprintf('%d:%d',hour,minute);
elseif hour < 10 && minute >=10
    time = sprintf('0%d:%d',hour,minute);
end
end

Github項目

完整源碼已上傳Github項目，點擊查看。

后記

這次是我第一次使用結構體來編程，感覺matlab真的很好用。我覺得作為建模編程人員一定要熟悉matlab中的矩陣、元胞數組、結構體、類的用法，能夠對數據按要求導入導出，並且模型代碼盡可能模塊化，減少冗余。