Python优化调度：基于改进遗传算法求解带时间窗约束多卫星任务规划_python遗传算法带有约束的问题

前言

参照网上基于改进遗传算法求解带时间窗约束多卫星任务规划方法，改进为基于Python实现并进行可视化的实现方法，实验的本质是一个带时间窗口的VRP的问题。

问题描述

假设全国现共有M个地面观测设备（每个观测设备都需要对卫星执行相应的观测任务），N个待观测卫星，且M<<N。每个待观测卫星相对于不同的地面设备都有P个可供选择的可见时间窗口。其中，每个可观测设备都可以在任何待观测卫星与之对应的可见时间窗口内对该卫星进行观测，观测的时长根据实际任务中该卫星所需观测时间而不同。同时，任意一颗卫星都可在其可见时间窗口内被地面观测设备所观测。由于观测设备自身的物理特性，每个地面观测设备对一颗卫星进行观测结束后对下一颗卫星观测之前，都需要经过设备转换时间，设备的转换时间根据设备的自身特性不同而不同。
简而言之，有M个地面观测设备，N个待观测卫星，需要为每个卫星指定地面站以及观测时间。

目标函数

总体观测时间尽量短

约束

1.每个卫星都需要被观测一次
2.地面站同时只能观测一个卫星
3.卫星需要在特定的时间窗口内，才能被地面站观测
4.卫星需要被观测足够长的时间

决策变量

1.卫星对应哪个地面站
2.卫星被观测的开始时间、结束时间

优化目标

卫星带时间窗的任务规划，类似于车间调度问题。对车间调度问题来说，存在多个工件，多台机床，每个工件又需要多道工序，每个工序只能在指定的机床上加工。

本文中的问题可以理解为车间调度问题的简化版，卫星对应工件，机床对应地面站。多道工序这里简化为卫星只需要被观测一次。同时该问题增加了一个时间窗，即卫星跟地面站只能在特定的时间段内可见，时间段长度、数量不定。

其实针对车间调度已经有很多成熟的编码方法，这里只介绍本文使用的编码方法。该编码分成两层，第一层是卫星的顺序编码，第二层是卫星对应的地面站编号。假设共有12个卫星，4个地面站，那么其中一个解可以表示如下：

图中第一个红框表示卫星的观测顺序，第二个框表示对应的地面站编号。只要确定了这两个，那么卫星对应的观测时间就可以计算出来。

基本思路为：卫星序列从头到尾遍历一遍，然后看当前卫星指定的地面站的时间窗口是什么，如果第一个时间窗口被别的卫星占用或者时间窗的持续时间不够观测卫星，那么就看看下一个时间窗满不满足约束。如果全部时间窗都不满足，那么计这个解违反约束，违反的程度为卫星需要被观测的时间。代码如下：

常见的使用遗传算法来处理实际工程问题，约束的处理可以使用罚函数方法。具体来说，将个体违反约束的程度，乘以一个系数加到目标函数中，那么违反约束的个体，会比满足约束的个体更加具有竞争力。从而推动算法向满足约束条件的方向进化。同时，遗传算法采用轮盘赌来选择下一代个体，具有更好的目标值的个体被选中的概率更大。

个体的违反约束情况作为第一排序准则，个体的适应度值作为第二排序准则。这种做法可以理解为可行性法则。具体来说，1）如果两个个体都违反约束，那么违反程度小的个体较好；2）如果两个个体一个违反、一个满足，那么满足约束的较好；3）如果两个个体都满足约束，那么，目标函数值小的个体较好。

在排序完成之后，删除种群中重复的个体（可能由变异产生）。这里可以理解为，决策变量完全一致的个体（在保持多样性上，还有其他方法可以选择，本文研究的问题是离散问题，判断冲不重复即可，其他情况可以判断个体之间的距离）。
挑选前N个个体组成下一代。

结果

代码

首先导入需要的库

import numpy as np
import pandas as pd
import os
import matplotlib.pyplot as plt
import matplotlib.patches as patches
1
2
3
4
5

定义全局变量

# 定义全局变量
Global = {}
Global['num_satellite'] = 4  # 地面站数量
Global['num_object'] = 12  # 观测目标数（卫星数）
Global['rank_satellite'] = None
Global['rank_object'] = None
Global['sat_need_time'] = None
Global['visible_window'] = {}
Global['num_visible_window'] = np.zeros((Global['num_object'], Global['num_satellite']))
# 读取数据
# 卫星优先级（列表/矩阵）
data_G = pd.read_csv('data/G.csv', encoding='ANSI').iloc[0:4, 1]
Global['rank_satellite'] = data_G.values
#print(Global['rank_satellite'])
# 观测目标优先级（列表/矩阵）
data_P = pd.read_csv('data/P.csv', encoding='ANSI').iloc[0:12, 1]
Global['rank_object'] = data_P.values
# 观测目标观测时长（列表/矩阵）
data_need = pd.read_csv('data/need.csv', encoding='ANSI').iloc[0:12, 1]
Global['sat_need_time'] = data_need.values

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

读取卫星数据

for i in range(0, Global['num_object']):
    datfile = f'data/sat{i + 1}.csv'
    data = pd.read_csv(datfile, header=None).iloc[0:4, 1:13].values
    for j in range(0, Global['num_satellite']):
        index = data[j] != 0
        Global['visible_window'][(i, j)] = data[j][index]
        # print( Global['visible_window'][(i, j)])
        Global['num_visible_window'][i, j] = len(Global['visible_window'][(i, j)]) // 2
        #  print(Global['num_visible_window'][i, j])
print('---test---')
1
2
3
4
5
6
7
8
9
10

初始化

def Init(N):
    empty = {'decs': None, 'objs': None, 'cons': None}
    population = [empty.copy() for _ in range(N)]

    for i in range(N):
        decs = np.concatenate([np.random.permutation(Global['num_object']) + 1,#卫星随机排列
                               np.random.randint(1, Global['num_satellite'] + 1, Global['num_object'])])
        # print(decs)
        population[i]['decs'] = decs.tolist()
        # print(population[i]['decs'])
    # Assuming there is a CalObj function to calculate objectives for the population
    population = CalObj(population)
    return population
1
2
3
4
5
6
7
8
9
10
11
12
13

def CalObj(population):
    N = len(population)
    for i in range(N):
        ind = population[i]
        # print(ind)
        object_list = ind['decs'][:Global['num_object']]
        # print(object_list)
        satellite_list = ind['decs'][Global['num_object']:]
        satellite_next_release_time = np.zeros(Global['num_satellite'])
        time_start_guance = np.zeros(Global['num_object'])
        # print(time_start_guance)
        time_end_guance = np.zeros(Global['num_object'])
        index_window_guance = np.zeros(Global['num_object'])
        cons = 0

        for j in range(Global['num_object']):
            # print((i,j))
            cur_object = object_list[j] - 1
            cur_satellite = satellite_list[j] - 1
            flag = 0
            for m in range(1, int(Global['num_visible_window'][cur_object, cur_satellite]) + 1):
                # if m==5:
                #     os.system("pause")
                time_start = Global['visible_window'][(cur_object, cur_satellite)][2 * m - 2]
                time_end = Global['visible_window'][(cur_object, cur_satellite)][2 * m - 1]
                # print((time_start, time_end))
                if satellite_next_release_time[cur_satellite] > time_end:
                    continue
                time_begin = max(satellite_next_release_time[cur_satellite], time_start)
                if time_begin < time_end - Global['sat_need_time'][cur_object]:
                    time_start_guance[cur_object] = time_begin
                    time_end_guance[cur_object] = time_start_guance[cur_object] + Global['sat_need_time'][cur_object]
                    satellite_next_release_time[cur_satellite] = time_end_guance[cur_object] + 60
                    index_window_guance[cur_object] = m
                    flag = 1
                    break
            if flag == 0:
                cons += Global['sat_need_time'][cur_object]
        T = max(time_end_guance)
        total_rank = 0
        for j in range(Global['num_object']):
            cur_object = object_list[j] - 1
            cur_satellite = satellite_list[j] - 1
            total_rank += Global['rank_satellite'][cur_satellite] * Global['rank_object'][cur_object]

        population[i]['objs'] = T - 10 * total_rank
        population[i]['cons'] = cons
        population[i]['time_start_guance'] = time_start_guance
        population[i]['time_end_guance'] = time_end_guance
        population[i]['satellite_list'] = satellite_list
        print(population[i]['satellite_list'])
        # 观测窗口编号
        population[i]['index_window_guance'] = index_window_guance
        population[i]['satellite_next_release_time'] = satellite_next_release_time
    return population
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55

交叉变异

def Mutate(population, state):
    # population = init(10)
    N = len(population)
    empty = {'decs': None, 'objs': None, 'cons': None}
    Offspring_temp = [empty.copy() for _ in range(N)]
    # Offspring = population.copy()
    for i in range(N):
        p1 = population[i]['decs'].copy()
        # print("交叉之前的offspring数据：\n",population[i]['decs'])
        if np.random.rand() < 0.8:  # 交叉
            # p2 = np.random.randint(N)
            while 1:
                p2_temp = np.random.randint(N)
                if p2_temp == i:
                    continue
                else:
                    break
                # p2 = np.random.randint(N)
            p2 = population[p2_temp]['decs'].copy()
            # 随机挑选两个进行片段交叉
            pos = np.sort(np.random.permutation(Global['num_object'])[:2])
            for j in range(pos[0], pos[1] + 1):
                if p1[j] != p2[j]:
                    ind = p1[:Global['num_object']].index(p2[j])
                    p1[ind] = p1[j]
                    p1[j] = p2[j]
            p1[pos[0] + Global['num_object']:pos[1] + Global['num_object'] + 1] = p2[pos[0] + Global['num_object']:pos[1] +Global['num_object'] + 1]

            if len(p1) < 24:
                print("测试交叉是否这里出现了问题：", len(p1), "第", i, "次")
                print("此时交叉的位置是:", pos)

        if np.random.rand() < 0.4:  # 变异
            pos = np.sort(np.random.permutation(Global['num_object'])[:2])  # 随机挑选两个位置进行片段逆转
            tmp = p1.copy()
            p1[pos[0]:pos[1] + 1] = np.flipud(tmp[pos[0]:pos[1] + 1])
            if len(p1) < 24:
                print("测试变异是否这里出现了问题：", len(p1), "第", i, "次")
                print("此时变异的位置是:", pos)
            pos = pos + Global['num_object']
            p1[pos[0]:pos[1] + 1] = np.flipud(tmp[pos[0]:pos[1] + 1])
        Offspring_temp[i]['decs'] = p1
    # for i_temp in range(N):
    #     Offspring[i_temp]['decs'][0:23] = Offspring_temp[i_temp]
    Offspring = CalObj(Offspring_temp)
    return Offspring
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46

挑选新个体

def Select(population, offspring, N):
    # 对每一代种群中的染色体进行选择，以进行后面的交叉和变异
    joint = population + offspring
    objs = np.array([ind['objs'] for ind in joint])
    cons = np.array([ind['cons'] for ind in joint])
    index = np.lexsort((objs, cons))
    joint = [joint[i] for i in index]
    # 删除重复个体
    del_indices = []
    for i in range(len(joint) - 1):
        if i in del_indices:
            continue
        for j in range(i + 1, len(joint)):
            if joint[i]['decs'] == joint[j]['decs']:
                del_indices.append(j)
    joint = [joint[i] for i in range(len(joint)) if i not in del_indices]
    population = joint[:N]
    return population
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

最后是主函数

if __name__ == '__main__':
    # 算法参数
    maxgen = 300
    popsize = 150
    population = Init(popsize)

    trace_obj = np.zeros(maxgen)
    trace_con = np.zeros(maxgen)
    # 进化开始
    for i in range(maxgen):
        # 交叉变异
        offspring = Mutate(population, i / maxgen)
        # 挑选新个体
        population = Select(population, offspring, popsize)
        # print(len(population))
        # 记录信息
        bestobj = population[0]['objs']
        trace_obj[i] = bestobj
        trace_con[i] = population[0]['cons']

        if not i % 10:
            cons = [ind['cons'] for ind in population]
            num = sum(1 for c in cons if c == 0)
            avgcons = np.mean(cons)
            print(f'第 {i} 代，满足约束个体数量：{num}，最佳个体：{bestobj}')

    # 进化结束
    print("---进化结束---\n---下面是绘图操作---")
    # 展示结果
    plt.plot(trace_obj) # 进化的迭代图
    # plt.title('最优目标值进化示意图')
    plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32