2012计算机仿真大赛（jsp 问题建模）

加工调度问题的遗传算法仿真

（matlab）

摘要

这是一个典型的Job-Shop动态排序问题。由于问题的连环嵌套性，使得用图解方法也变得不切实际。传统的运筹学方法，即便在单目标优化的静态调度问题中也难以有效应用。

本文独辟蹊径，开创性地提出了加工调度问题的时间延时器模型，在此基础上采用现代进化算法中有代表性发展优势的遗传算法求解，文章有针对性地选取遗传算法关键环节的适宜方法，采用MATLAB软件实现算法模拟，得出优化方案，并与计算机随机模拟结果加以比较显示出遗传算法之优化效果。

首先，我们考虑到最简单的情况，在T=2×2（两个工件，两台机器），T=3×3（三个工件，三台机器），情况下，给出仿真结果，并用常规算法得出最优解，验证遗传算法的优化性质。

T=2×2（2个工件，2台机器）； T=2×3（2个工件，3台机器）；

在此基础上我们增大T的维数进行高维仿真：

T=10×10（10个工件，10台机器） T=100×10（100个工件，10台机器）

维数T 总加工时间t T=2×2 158 T=2×3 96 T=10×10 550 T=100×10 850 文章有针对性地选取遗传算法关键环节的适宜方法，采用MATLAB软件实现算法模拟，得出优化方案，并与计算机随机模拟结果加以比较显示出遗传算法之优化效果。 对车间调度系列问题的有效解决具有一定参考和借鉴价值。

一．问题重述

工厂中，有n个不同的配件需要生产，每个配件都必须由m台不同的机器进行顺序加工处理，配件i在机器j上所需的处理时间为t(i,j)。现约定未完工前不允许中断处理，配件不能拆分成更小配件。

条件：1、每件产品必须按规定的工序加工，不得颠倒； 2、每台设备在同一时间只能担任一项任务。 3、配件不能拆分成更小配件。

问题：做出生产安排，希望在尽可能短的时间里，完成所接受的全部任务。 要求：给出每台设备承担任务的时间表。

二．模型假设

1．每一时刻，每台机器只能加工一个工件，且每个工件只能被一台机器所加工 ，同时加工过程为不间断； 2．所有机器均同时开工，且工件从机器I 到机器J 的转移过程时间损耗不计； 3．各工件必须按工艺路线以指定的次序在机器上加工多次； 4．操作允许等待，即前一操作未完成，则后面的操作需要等待，可用资源有限。

三．符号说明及初始数据表达分析

Ji - 第i个工件 （i=1?6）

（i,j）JM- 机器顺序阵 JM表示i工件的第 j个操作的机器号

Mj- 第j台机器 （j=1?4）

MJ- 工件排列阵 MJ(i,j)表i机器上第j次加工的工件号 T - 加工时间阵 T(i,j)为i工件的第 j个操作的时间周期 C - 整个任务完成时间 H - 延时器 Ht -延迟时间

四．模型提出

延时器（H）：延时器是通信与信息系统中提出的一种理想化模型，用来描 述信号通过一个系统时的延时效应；

- 1 -

延迟时间（Ht）：信号通过延时器时，在时间域上的延迟时间；

对于一个用N个代加工产品的系统，假设有M台机器，则每台机器可以等价为一个延时器H，且延时时间 Ht（i，j），那么每一个产品i在机器j上的加工过程可以建模成信号通过系统Ht（i，j）的过程。下面通过这一理论进行研究。

如果按照一定规则排序，当多个工件出现“抢占”同一机器的局面的时候,我们可以制定如下的工序安排规则:

1. 优先选择总剩余时间或总剩余操作较多的工件。（如果出现总剩余加工时间多者总剩余操作数反而较少的情况时,按照程度具体情况具体分析）。 2. 机器方面来说,尽量避免等待空闲时间，优先考虑剩余净加工时间或者剩余加工总次数较多的机器,尤其是机器 D ,即倘若能够使机器D不间断工作且其他机器完工时间均不多余75时,那么就可以得到最优解 。

首先按照最优化时间为75的设想避免D出现等待，排序后得到升以下具体排列顺序。

各机器承担任务表为(其中粗体字为对应工件产品号，括号内为对应时间周期段)： 操作1 操作 2 操作3 操作4 操作5 操作6 操作7 操作8 操作9 操作10 6 2 1 6 3 4 5 6 3 6 (1) (2-5) (12-13) (14-20) (21-35) (36) (43-54) (55-56) (57-64) (66-73) A B C D 4 (1-7) 6 1 3 5 4 4 1 5 5 3 6 6 2 1 2 (8-11) (12-15) (16-19) (36-55) (56-58) 3 (1-3) 5 4 6 1 2 5 (75) 4 (4-11) (12-17) (18-25) (26-28) (29-52) (55-60) (61-65) (66-69) (70-72) 5 (表四) A机器14B机器C机器D机器03101078^424672074815^321(1-10) (11-17) (18-38) (39-42) (43-47) (48-52) (53-68) (69-74) 1^2024451223^68^548367018055016603040 - 2 -

(图一)

上图为加工周期图（甘特图），标注数字为相应操作的周期，完工时间为第75周期。

五．计算机随机模拟（编程）

1．编码：

随机产生生产的工序操作优先顺序，进行编码，如：K=[ 4 3 5 6 6 2 3 1 4 1 6 3 5 4 5 3 6 6 4 1 5 5 1 3 2 6 2 2 4 4 1 5 6 6 5 ] （注：同时作为下文的染色体之用） 意思为：工件4优先被考虑进行第一次操作，然后3进行其第一步操作，然后5操作，6操作，再6操作其第二步工序，依次进行。如果前后互相不冲突，则可同时在不同机器上操作。

通过排列组合得出，总共有类似K的排列序列 2?1023多种！

当然，这其中只对应解 [75，247]，意味着有大量排列序列对应同一加工方案，而大量加工方案又对应同一时间解。

2．解码：

即对编码进行翻译，产生具体可操作工序安排方案，这里采用活动化

解码算法。例如工件2第i步操作（记为JM（2，i），且在机器A上进行）被安排在工件3第j步操作（记为JM（3,j））后面进行，那么如果安排好JM （3,j）后，只要JM（2，i）在工件2已经排序好的操作之后进行，那么操 作JM（2，i）可插入到机器A处最前可安置的时间段进行。

在这里，一个编码序列对应一个加工方案，而一个加工方案可对应一个或多个编码序列，这就是二者之关系。

3．编程：

通过一组随机编码产生一生产加工优先序列，通过解码过程产生相应加

工方案及其总耗费时间C . N次模拟后即可得出解C的概率密度分布情况以及相对最优解（N个C的最小值，如80，77等，甚至出现75）。

4．计算机模拟所得数据分析

a. 进行100次模拟得出最优解情况： (共运行10次)

82，83，82，84，78，80，81，83，87，82 平均值 82.2，每回耗时约3秒 b. 进行1000次模拟得出最优解情况：(共运行10次)

80，79，78，78，79，79，76，80，77，78 平均值 78.4 , 每回耗时25秒 c. 进行10000次模拟得出最优解情况：(共运行10次)

76，77，77，75，76，76，77，76，76，77 平均值 76.3, 每回耗时4分钟

- 3 -

上( 图二 ) 算法的复杂度。

结论：

从图二中可以看出，由于问题的连环嵌套性，使得用图解方 法也变得不切实际。

23 如果想将2?10中排序序列以平均出现一次的可能性进行模拟，

23即运行2?10次，计算机需运行将近150万亿年！当然，我们没有这个必要，因为我们只需要运行数万次，就很可能得到最优解75。

六．遗传算法模型建立和步骤解法

遗传算法（Genetic Algorithm）作为一种优化算法特别适合于对象模型难于建立、搜索空间非常庞大的复杂问题的优化求解。它和模糊控制技术一样，虽然在理论上还没有完善，但是在实践中已经得到了广泛的应用。遗传算法的基本思想是：模仿生物系统“适者生成\的原理，通过选择、复制、交叉、变异等简单操作的多次重复来达到去劣存优的目的，从而获得问题的优化结果。遗传算法的实现由两个部分组成，一是编码与解码，二是遗传操作。其中遗传操作又包括选择、复制、交叉、变异等步骤。

本文根据实际情况采取了1-6整数编码。数字1，2，3，4，5，6分别代表6件待加工产品。

本文遗传算法基本流程:

1、初始化参数： 族群 60，循环500次，交叉0.8 变异 0.6 代沟0.9

2、初始化群：按调度优先级编码，比如 3个零件，每个零件3个工序，就初始化：

1、3、4、5、6、7、8、9、2 2、1、3、4、5、6、7、8、9 等等

3、计算适应值：解码成工序序列，计算时间

4、选择：在原族群中，按轮盘法选择出 60×0.9（代沟）=54个个体组成 新族群。

- 4 -

5、交叉： 对选择出来的新族群进行交叉，在新族群中，随机选择2个个体（选择过的个体将不在选择，保证所有个体都选择），随机生成有一个数，如果大于交叉概率，就进行2点交叉。2点交叉位置每次都是随机的

比如 1、2、3、5、6、7、8、4、9 2、1、3、5、6、4、9、7、8 如果交叉位置 在 2和5

先 生成 0、2、3、5、6、0、0、0

然后按 下面的数据 删除2、3、5、6后写入 就是 1、2、3、5、6、4、9、7、8

6、变异：对交叉出来的新族群进行变异，在新族群中，对每个个体，随机生成有一个数，如果大于变异概率，随机生成2个位置数据，将2个位置的数据进行交换

7、替换：新的族群个体是54个，原来的族群为60，保留6个适应值好的，其他代替54个。

三、计算效果图

3500解的变化种群均值的变化 300025002000150010005000 050100150200250300350400450500

1098765432102090706080827150408172083521039421313324400600611048896105478001510826161000683866991200140010010174459254852332951146932276628777346463351213252757986759697975444386106891436107487828376597538410949175618193941301024291735120029585510

甘特图

通过编码，解码程序随机产生N个（有一定数量,如50或100）个体构成初

始种群

a) 从初始中群中选取2个具有最优染色体（最有排序方案）的个体作为临时个

体（父代）；

b) 如果此2个体中有一个个体通过解码操作能够实现最优排序（即使总时间为

75周期），那么结束此算法，得到最优解；

c) 对2个临时个体以一定方式(循环交叉)执行染色体交叉变换和变异选择（小

概率,互换操作），产生2个新的个体；

- 5 -

d) 对父代和子代共4个个体进行选择，从中选出最佳的2个个体，做为下一代

的父代；

e) 重复执行第二步(b)操作；

f) 如果执行完M步后仍然未得出答案75，那么将目前的最优解作为本算法的最

优解答案。 1．编码和解码 （同上）

2．交叉变换(crossover)

对2个父代临时个体进行染色体交叉变换，采用循环交叉方法（Cycle crossover CX），如父代染色体为：X：[9 2 6 4 7 3 5 8 1]和Y：[3 4 5 8 1 6 7 2 9]，如果随机选到第二位开始交叉，那么X的2对应Y的4，X的4对应Y的8，X的8对应Y的2，这样就确定了以上为不变的染色体，其余位置的染色体互换位置，最后得到X?: [3 2 5 4 1 6 7 8 9], Y?: [9 4 6 8 7 3 5 2 1],实现交叉变换。

3．变异选择（mutation）

采用互换操作（SWAP），，即随机交换染色体中两不同基因的位置。如上面的染色体为：X?: [3 2 5 4 1 6 7 8 9] 。随机产生变换位置号，如产生随机数3和5，那么交换数字后得到染色体： [3 2 1 4 5 6 7 8 9]， 变异概率取0.1 。

4．选择操作（selection）

对父代2个体f1,f2和子代2个体f3,f4进行选择，通过编码操作确定

具有最优解的2个个体，成为新一代f1和f2 。

如此，通过多次编码和解码随机产生一定数量的个体，选取2个最佳个体进行交叉变换操作，产生2个新个体，然后对4个个体进行选择，产生下一代，如果某时刻通过解码操作得出最优解（所有解的下限，这里是75周期），那么算法结束，否则循环进行，直至预先给定的循环次数达到为止，以最后得到的最优解作为最终最优解。

七．遗传算法仿真（采用MATLAB工具编程）

主程序如下：（子程序见略）

clc clear

close all

%调度时间矩阵

%产生机器生产周期矩阵t(i,j),第i件产品在地j台机器上的生产时间 %task=('产生生产周期矩阵t（i，j）');

%T=input('输入加工周期矩阵t(i,j),第i件产品在地j台机器上的生产时间'); %[m,n]=size(T)

- 6 -

a=input('输入要加工产品个数'); b=input('输入工序数目'); T=randint(a,b,[0,100]) [m,n]=size(T)

M=500;%input('输入遗传进化迭代次数：');遗传进化迭代次数 N=30;%input('输入种群规模，取偶数：');种群规模，偶数 Pm=0.02;%input('输入变异概率：')变异概率

P=ones(1,n);%[1 1 1 1 1 1 1 1 1 1];% m道工序加工时每到工序都由一台机器加工；

[Zp,Y1p,Y2p,Y3p,Xp,LC1,LC2]=gajsp(M,N,Pm,T,P) % JSPGA.m

% 车间作业调度问题遗传算法 % 输入参数列表

% M 遗传进化迭代次数 % N 种群规模(取偶数) % Pm 变异概率

% T m×n的矩阵，存储m个工件n个工序的加工时间

% P 1×n的向量，n个工序中，每一个工序所具有的机床数目 % 输出参数列表

% Zp 最优的Makespan值

% Y1p 最优方案中，各工件各工序的开始时刻，可根据它绘出甘特图 % Y2p 最优方案中，各工件各工序的结束时刻，可根据它绘出甘特图 % Y3p 最优方案中，各工件各工序使用的机器编号

% Xp 最优决策变量的值，决策变量是一个实数编码的m×n矩阵

% 最后，程序还将绘出三副图片：两条收敛曲线图和甘特图（各工件的调度时序图）

八．遗传算法仿真结果

首先给出最优方案:

在进行某次n=2,m=2的操作后得到仿真最优结果75周期时间：

T =[22 56； 48 75]

各工件起始加工时间： Y1p =

0 63 63 163

各工件结束时时间： Y2p =

63 126 163 163

- 7 -

决策变量值： Xp =

1.4253 1.9866 1.3377 1.4705

minminmax = 163 （最终最优解）

进行n=5,m=4的操作后得到仿真最优结果75周期时间

T =

98 61 4 75 66 95 74 61 14 51 33 53 10 28 54 75 75 25 4 59 各工件起始加工时间： Y1p =

0 98 159 163 238 98 193 193 238 193 193 226 226 226 244 226 301 301 301 301 各工件结束时时间： Y2p =

98 159 163 238 238 193 193 193 193 244 226 226 226 226 226 301 301 301 301 301 决策变量值： Xp =

1.6832 1.8196 1.5654 1.2095 1.0617 1.3252 1.5940 1.0633 1.6320 1.8301

- 8 -

1.8724 1.9839 1.3120 1.7488 1.3415 1.0092 1.2842 1.0597 1.5039 1.1822 minminmax = 301 （最终最优解）

首先给出最优方案:

在进行某次n=10,m=10的操作后得到模拟最优结果75周期时间：

T =

86 97 7 24 90 70 68 37 12 0 0 80 22 25 49 64 48 90 10 79 75 17 2 45 94 26 31 48 62 33 10 72 66 58 96 79 86 28 7 97 18 94 12 26 21 20 49 31 67 31 70 63 8 67 71 89 4 30 88 40 65 5 24 42 5 3 84 42 34 90 74 66 41 76 87 4 19 26 63 94 56 76 14 64 36 28 66 5 46 73 79 17 95 10 70 3 14 32 57 11

各工件起始加工时间： Y1p =

0 86 183 183 183 183 183 272 183 195 86 183 86 108 133 182 257 182 272 272 86 161 161 161 182 161 187 187 195 187 161 171 171 171 171 187 171 257 257 257 171 189 189 189 189 189 189 189 189 220 189 259 259 259 259 259 259 259 259 259 259 324 324 324 324 324 324 324 324 324 324 398 398 398 398 398 398 398 398 398 398 454 454 454 454 454 454 454 454 454 454 533 533 533 533 533 533 533 533 533

- 9 -

各工件结束时时间： Y2p =

86 183 183 183 183 183 183 183 195 195 86 86 108 133 182 182 182 272 272 272 161 161 161 161 161 187 187 187 187 220 171 171 171 171 171 171 257 257 257 257 189 189 189 189 189 189 189 189 189 189 259 259 259 259 259 259 259 259 259 259 324 324 324 324 324 324 324 324 324 324 398 398 398 398 398 398 398 398 398 398 454 454 454 454 454 454 454 454 454 454 533 533 533 533 533 533 533 533 533 533

决策变量值： Xp =

1.0513 1.4660 1.2732 1.3301 1.7228 1.3805 1.3043 1.7735

1.9758 1.0491 1.9426 1.9652 1.7827 1.6814 1.3187 1.2158

1.8303 1.5809 1.8688 1.8334 1.7009 1.6825 1.1696 1.9176

1.6436 1.6149 1.8296 1.6525 1.5913 1.9028 1.1321 1.3234

1.7964 1.3149 1.5326 1.8349 1.0311 1.9893 1.7632 1.8012

1.5285 1.4037 1.9852 1.5943 1.2874 1.9002 1.3755 1.6691

1.1942 1.2464 1.1947 1.7714 1.1543 1.4585 1.5075 1.2441

1.6136 1.9751 1.0957 1.7946 1.8589 1.0098 1.2426 1.7188

1.7506 1.6793 1.6136 1.9287 1.3694 1.3961 1.9027 1.7613

1.5813 1.4417 1.0425 1.4156 1.6768 1.0500 1.9479 1.9695

minminmax = 533 （最终最优解）

- 10 -

1.4054 1.6063 1.9691 1.3525 1.0028 1.6685 1.1436 1.3321 1.0600 1.8223 1.9396 1.2957 1.0763 1.2238 1.6558 1.0802 1.3175 1.0953 1.1725 1.9244

以下为:取不同N和M值情况下数据优化过程以及时间上的比较: ( 表五 ) N 2 2 2 10 10 100 100 M 2 10 100 10 100 100 1000 第一次运行 第二次运行 第三次运行 第四次运行 第五次运行 平均运行时间 120 256 825 922 1025 1542 2024 131 284 822 929 1020 1548 2048 122 250 826 925 1010 1550 2052 132 280 828 930 1023 1561 2047 125 266 830 916 1035 1547 2096 / / / / / / 11(s) 明显地发现: M 的增大所产生的优化效果明显好于N增大产生的优化效果

M 从1-10-100-1000-10000的变化使最优解有从9*--8*--7* 的变化规律， N 的1-10-100的变化则不明显。

因此相对于随机取解, 经过遗传算法优化之后效果是显著的。

另外对采用遗传算法前后模拟所得数据进行比较: 第一次 第二次 第三次 第四次 第五次 均值 平均时间 N=1000 ,M=0 78 79 81 80 79 79．4 25(S) N=1,N=1,M=1000 80 75 77 76 82 78.0 8.5(S) (表六) 由上表看来,虽然在均值方面相差不显著,但是时间上采用遗传算法之后节约了约三分之二的运行时间,效率显而易见。

九．模型优缺点及改进

模型优点，采用遗传算法可对一个理论上无法求尽的NP问题以较有效方法在较短时间内得到较精确解。

缺点，采用遗传算法还不全面，只是一个简单模型，未考虑收敛性,适配

- 11 -

值等，对参数设置与最优解关系研究还不够深入。 模型本身还具有漏动，仍需改进，较好设想为采用和模拟退火算法的混合遗传算法，由于时间紧迫，在此就不再深入给出模型及分析了。

[参考文献]：

[1].车间调度与遗传算法 王凌 清华大学出版社 [2].数值计算的算法与分析 张可村 赵英良 科学出版社

[3].Permutation Based GAs and Ordered Greed Peter G. Anderson, [4].MATLAB6.0 与科学计算 王沫然 电子工业出版社 [5].C程序设计（第二版） 潭浩强 清华大学出版社

附录：

主程序：

clc clear

close all

%调度时间矩阵

%产生机器生产周期矩阵t(i,j),第i件产品在地j台机器上的生产时间 %task=('产生生产周期矩阵t（i，j）');

%T=input('输入加工周期矩阵t(i,j),第i件产品在地j台机器上的生产时间'); %[m,n]=size(T)

a=input('输入要加工产品个数'); b=input('输入工序数目'); T=randint(a,b,[0,100]) [m,n]=size(T)

%T=[29 78 9 36 49 11 62 56 44 21; % 43 90 75 11 69 28 46 46 72 30; % 91 85 39 74 90 10 12 89 45 33; % 81 95 71 99 9 52 85 98 22 43; % 14 6 22 61 26 69 21 49 72 53; % 84 2 52 95 48 72 47 65 6 25; % 46 37 61 13 32 21 32 89 30 55;

- 12 -

% 31 86 46 74 32 88 19 48 36 79; % 76 69 76 51 85 11 40 89 26 74; % 85 13 61 7 64 76 47 52 90 45];

%Jm=[0 1 2 3 4 5 6 7 8 9; % 0 2 4 9 3 1 6 5 7 8; % 1 0 3 2 8 5 7 6 9 4; % 1 2 0 4 6 8 7 3 9 5; % 2 0 1 5 3 4 8 7 9 6; % 2 1 5 3 8 9 0 6 4 7; % 1 0 3 2 6 5 9 8 7 4; % 2 0 1 5 4 6 8 9 7 3; % 0 1 3 5 2 9 6 7 4 8; % 1 0 2 6 8 9 5 3 4 7];

M=50;%input('输入遗传进化迭代次数：');遗传进化迭代次数 N=30;%input('输入种群规模，取偶数：');种群规模，偶数 Pm=0.02;%input('输入变异概率：')变异概率

P=ones(1,n);%[1 1 1 1 1 1 1 1 1 1];% m道工序加工时每到工序都由一台机器加工；

[Zp,Y1p,Y2p,Y3p,Xp,LC1,LC2]=gajsp(M,N,Pm,T,P) T

%-------------------------------------------------------------------------- % JSPGA.m

% 车间作业调度问题遗传算法

%--------------------------------------------------------------------------

% 输入参数列表

% M 遗传进化迭代次数 % N 种群规模(取偶数) % Pm 变异概率

% T m×n的矩阵，存储m个工件n个工序的加工时间

% P 1×n的向量，n个工序中，每一个工序所具有的机床数目 % 输出参数列表

% Zp 最优的Makespan值

% Y1p 最优方案中，各工件各工序的开始时刻，可根据它绘出甘特图 % Y2p 最优方案中，各工件各工序的结束时刻，可根据它绘出甘特图 % Y3p 最优方案中，各工件各工序使用的机器编号

% Xp 最优决策变量的值，决策变量是一个实数编码的m×n矩阵 % LC1 收敛曲线1，各代最优个体适应值的记录 % LC2 收敛曲线2，各代群体平均适应值的记录

% 最后，程序还将绘出三副图片：两条收敛曲线图和甘特图（各工件的调度时序图）

- 13 -

程序流程图： 子程序一：

function [Zp,Y1p,Y2p,Y3p]=COST(X,T,P,plotif) % JSPGA的内联子函数，用于求调度方案的Makespan值 % 输入参数列表

% X 调度方案的编码矩阵，是一个实数编码的m×n矩阵 % T m×n的矩阵，存储m个工件n个工序的加工时间

% P 1×n的向量，n个工序中，每一个工序所具有的机床数目 % plotif 是否绘甘特图的控制参数 % 输出参数列表

% Zp 最优的Makespan值

% Y1p 最优方案中，各工件各工序的开始时刻 % Y2p 最优方案中，各工件各工序的结束时刻 % Y3p 最优方案中，各工件各工序使用的机器编号

%第一步：变量初始化 [m,n]=size(X); Y1p=zeros(m,n); Y2p=zeros(m,n); Y3p=zeros(m,n);

%第二步：计算第一道工序的安排 Q1=zeros(m,1); Q2=zeros(m,1);

R=X(:,1);%取出第一道工序

Q3=floor(R);%向下取整即得到各工件在第一道工序使用的机器的编号 %下面计算各工件第一道工序的开始时刻和结束时刻 for i=1(1)%取出机器编号

pos=find(Q3==i);%取出使用编号为i的机器为其加工的工件的编号 lenpos=length(pos); if lenpos>=1 Q1(pos(1))=0;

Q2(pos(1))=T(pos(1),1); if lenpos>=2

for j=2:lenpos

Q1(pos(j))=Q2(pos(j-1));

Q2(pos(j))=Q2(pos(j-1))+T(pos(j),1); end end end end

Y1p(:,1)=Q1; Y2p(:,1)=Q2;

- 14 -

本文来源：https://www.bwwdw.com/article/esnh.html