Matlab代码及解析（可运行）基于改进遗传算法求解带时间窗约束多卫星任务规划_带约束遗传算法matlab

作者：IT小白 | 2024-07-08 09:30:13

踩

带约束遗传算法matlab

问题介绍

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

简而言之，有M个地面观测设备，N个待观测卫星，需要为每个卫星指定地面站以及观测时间。

目标函数 是使得观测结束的时间最短
约束条件 可以表示为：

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

具体数学模型较为简单，不做展示。
决策变量： 1）卫星对应哪个地面站，2）卫星被观测的开始时间、结束时间

问题难点

使用遗传算法对该问题进行优化求解，存在的难点主要为：

如何对该问题进行编码（怎么表示一个解）
怎么处理问题中存在的约束条件
如何确保收敛性。一方面要保持种群的多样性，防止出现聚堆现象导致进化停滞不前；另一方面，如何设计合理的交叉变异方法。

以上问题基本就是使用进化算法求解实际工程问题中存在的共性问题。

分析

编码方法介绍

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

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

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

对应的，遗传算法种群编码函数为

// Init.m
function population = Init(N)

global Global
empty.decs = [];
empty.objs = [];
empty.cons = [];

population = repmat(empty,1,N);
for i=1:N
    population(i).decs = [randperm(Global.num_satellite,Global.num_satellite) ...
     randi([1,Global.num_ground],1,Global.num_satellite)];
end
population = CalObj(population);
1
2
3
4
5
6
7
8
9
10
11
12
13
14

*randperm(Global.num_satellite,Global.num_satellite)*表示产生卫星顺序的序列，*randi([1,Global.num_ground],1,Global.num_satellite)*表示指定的地面站编号

种群的数据结构表示为
在这里插入图片描述

目标函数的计算

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

// CalObj.m
function population = CalObj(population)
% population = Init(1);

global Global

N = length(population);

for i=1:N
    ind = population(i);
    satellite_list = ind.decs(1:Global.num_satellite);
    ground_list = ind.decs(Global.num_satellite+1:end);
    
    ground_next_release_time = zeros(1,Global.num_ground); %地面站的下一次可观测时间
    
    time_start_guance = zeros(1,Global.num_satellite); %开始观测时间
    time_end_guance = zeros(1,Global.num_satellite); %结束观测时间
    index_window_guance = zeros(1,Global.num_satellite); %观测窗口的编号
    cons = 0; %违反约束的情况
    
    for j = 1:Global.num_satellite
        cur_satellite = satellite_list(j); %当前卫星
        cur_ground = ground_list(j); %当前卫星选择的地面站
        
        % 依次检查当前卫星的观测时间窗
        flag = 0;
        for m = 1:Global.num_visible_window(cur_satellite,cur_ground)
            time_start = Global.visible_window{cur_satellite,cur_ground}(2*m-1); %观测窗口的开始时间
            time_end = Global.visible_window{cur_satellite,cur_ground}(2*m); %观测窗口的结束时间
            if ground_next_release_time(cur_ground) > time_end
                continue;
            end
            time_begin = max(ground_next_release_time(cur_ground),time_start);
            
            if time_begin < time_end - Global.sat_need_time(cur_satellite) %当前观测窗口可以用来观测
                time_start_guance(cur_satellite) = time_begin;
                time_end_guance(cur_satellite) = time_start_guance(cur_satellite) + Global.sat_need_time(cur_satellite);
                ground_next_release_time(cur_ground) = time_end_guance(cur_satellite) + 60; %设备的反应时间
                index_window_guance(cur_satellite) = m;
                flag=1;
                break;
            end
        end
        if flag == 0 %说明这个卫星没有被安排到地面站
            cons = cons + Global.sat_need_time(cur_satellite);
        end
    end
    
    %% 计算目标函数
    T = max(time_end_guance);
    total_rank = 0;
    for j = 1:Global.num_satellite
        cur_satellite = satellite_list(j); %当前卫星
        cur_ground = ground_list(j); %当前卫星选择的地面站
        total_rank = total_rank + Global.rank_ground(cur_ground) * Global.rank_satellite(cur_satellite);
    end
    
    %% 传出值
    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).ground_list = ground_list;
    population(i).index_window_guance = index_window_guance;
    population(i).ground_next_release_time = ground_next_release_time;
end
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
56
57
58
59
60
61
62
63
64
65
66

个体变异方法介绍

采用上面所说的编码方式，那么第一层编码非常类似常见的旅行商问题，编码表示为一串序列。那么类比于TSP问题，本文设计了两种变异方法。

随机挑选另一个个体，然后随机选择两个位置，将另一个个体两个位置之间的编码付给当前个体的相同位置；
随机挑选两个位置，将两个位置之间的编码逆转。

具体代码表示为：

// Mutate.m
function offspring=Mutate(population,state)
% 本函数完成变异操作
% population = Init(10);

global Global
N=length(population);

offspring = population;
for i=1:N
    p1 = population(i).decs;
    if rand < 0.8 % 交叉
        p2 = i;
        while p2==i
           p2 = randperm(N,1); 
        end
        p2 = population(p2).decs;
        pos = sort(randperm(Global.num_satellite,2)); %随机挑选两个位置进行片段交叉
        for j=pos(1):pos(2)
           if p1(j) ~= p2(j)
               ind = find(p1(1:Global.num_satellite)==p2(j));
               p1(ind) = p1(j);
               p1(j) = p2(j);
           end
        end
        
        p1(pos(1)+Global.num_satellite:pos(2)+Global.num_satellite) = p2(pos(1)+Global.num_satellite:pos(2)+Global.num_satellite);
    end
    
    if rand < 0.4 % 变异
        pos = sort(randperm(Global.num_satellite,2)); %随机挑选两个位置进行片段逆转
        tmp = p1;
        p1(pos(1):pos(2)) = tmp(pos(2):-1:pos(1));
        pos = pos + Global.num_satellite;
        p1(pos(1):pos(2)) = tmp(pos(2):-1:pos(1));
    end
    
    offspring(i).decs = p1;
end
offspring = CalObj(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

约束处理方法、种群收敛性及多样性

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

经典的算法存在很多问题，例如

罚函数的系数怎么设置？如果约束很难被满足，那么算法一直处于在不满足约束的情况下优化，得到的解没有用。
轮盘赌的选择策略，导致每一代最优秀的个体可能没被选上；或者经过很多次的迭代，发现种群中的个体全部一样，那么丧失了多样性，算法就停滞不前。对于第一个问题，有精英保留策略，虽有一定效果，但是保存多少个精英个体合适呢？对于第二个问题，轮盘赌之前的适应度值怎么计算是一个问题。

对于这个问题，改进的遗传算法作出如下调整：

个体的违反约束情况作为第一排序准则，个体的适应度值作为第二排序准则。这种做法可以理解为可行性法则。具体来说，1）如果两个个体都违反约束，那么违反程度小的个体较好；2）如果两个个体一个违反、一个满足，那么满足约束的较好；3）如果两个个体都满足约束，那么，目标函数值小的个体较好。
在排序完成之后，删除种群中重复的个体（可能由变异产生）。这里可以理解为，决策变量完全一致的个体（在保持多样性上，还有其他方法可以选择，本文研究的问题是离散问题，判断冲不重复即可，其他情况可以判断个体之间的距离）。
挑选前N个个体组成下一代。

具体代码如下：

// Select.m
function population=Select(population,offspring,N)
% 本函数对每一代种群中的染色体进行选择，以进行后面的交叉和变异

joint = [population,offspring];
objs = [joint.objs]';
cons = [joint.cons]';

[~,index] = sortrows([cons,objs]);

joint = joint(index);

% 删除重复个体
del = [];
for i=1:length(joint)-1
    if find(i==del)
        continue;
    end
   for j=i+1:length(joint)
       if isequal(joint(i).decs,joint(j).decs)
          del = [del j]; 
       end
   end
end
joint(del) = [];

population = joint(1:N);
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

程序主体

采用基本的遗传算法流程，代码如下：

// RunMe.m
%% 清空环境
clc
clear
close all

%% 模型参数
global Global
Global.num_ground = 4;    %设备数
Global.num_satellite = 12;    %卫星数

%% 读取数据
[a,~,~] = xlsread('data\G.csv');
Global.rank_ground = a';%设备优先级（列表/矩阵）
[a,~,~] = xlsread('data\P.csv');
Global.rank_satellite = a';%卫星优先级（列表/矩阵）
[a,~,~] = xlsread('data\need.csv');
Global.sat_need_time = a'*1.2;%卫星观测时长（列表/矩阵）

Global.visible_window = cell(Global.num_satellite,Global.num_ground);
Global.num_visible_window = zeros(Global.num_satellite,Global.num_ground);
for i=1:Global.num_satellite
    datfile = ['data\sat' num2str(i) '.csv'];
    [a,~,~] = xlsread(datfile);
    
    for j=1:Global.num_ground
        index = a(j,:)~=0;
        Global.visible_window{i,j} = a(j,index);
        Global.num_visible_window(i,j) = numel(Global.visible_window{i,j})/2;
    end
end

%% 算法参数
maxgen = 1;
popsize = 150;
population = Init(popsize);

trace_obj = zeros(1,maxgen);
trace_con = zeros(1,maxgen);

%% 进化开始
for i=1:maxgen
    % 交叉变异
    offspring = Mutate(population,i/maxgen);
    % 挑选新个体
    population = Select(population,offspring,popsize);
    
    % 记录信息
    bestobj = population(1).objs;
    trace_obj(i) = bestobj;
    trace_con(i) = population(1).cons;
    
    if ~mod(i,10)
        cons = [population.cons];
        num = sum(cons==0);
        avgcons = mean(cons);
        disp(['第' num2str(i) '代，满足约束个体数量：' num2str(num), '，最佳个体：' num2str(bestobj)])
    end
end
%进化结束

%% 展示结果
figure
plot(trace_obj)
title('最优目标值进化示意图')

bestsol = population(1);
drawresult
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
56
57
58
59
60
61
62
63
64
65
66
67
68

结果展示

代码如下：

// drawresult.m
close all

figure
global Global

% c_space = linspecer(Global.num_satellite);
c_space = colormap(jet(12));
for i=1:Global.num_satellite
    cur_satellite = bestsol.decs(i);
    cur_ground = bestsol.ground_list(i);
    ind_window = bestsol.index_window_guance(cur_satellite);
    
    subplot(4,3,cur_satellite)
    
    t_s = bestsol.time_start_guance(cur_satellite);
    t_e = bestsol.time_end_guance(cur_satellite);
    
    if t_s == 0 && t_e ==0
        continue;
    end
    
    t_s_window = Global.visible_window{cur_satellite,cur_ground}(2*bestsol.index_window_guance(cur_satellite)-1);
    t_e_window = Global.visible_window{cur_satellite,cur_ground}(2*bestsol.index_window_guance(cur_satellite));
    rec = [t_s_window,cur_ground-0.1,t_e_window-t_s_window,0.2];
    rectangle('Position',rec,'LineWidth',0.5,'LineStyle','-','FaceColor',c_space(i,:),'Curvature',0.5);%draw every rectangle
%     for j=1:Global.num_visible_window(cur_satellite,cur_ground)
%         t_s_window = Global.visible_window{cur_satellite,cur_ground}(2*j-1);
%         t_e_window = Global.visible_window{cur_satellite,cur_ground}(2*j);
%         rec = [t_s_window,cur_ground-0.1,t_e_window-t_s_window,0.2];
%         rectangle('Position',rec,'LineWidth',0.5,'LineStyle','-','FaceColor',c_space(i,:),'Curvature',0.5);%draw every rectangle
%     end
    
    rec = [t_s,cur_ground-0.25,t_e-t_s,0.5];
    rectangle('Position',rec,'LineWidth',0.5,'LineStyle','-','FaceColor',c_space(i,:),'Curvature',0.5);%draw every rectangle
    text(t_s+100,cur_ground,num2str(cur_satellite),'FontWeight','Bold','fontsize',8);
    ylim([0,5])
    title(['Satellite' num2str(cur_satellite)])
end

figure
for i=1:Global.num_satellite
    cur_satellite = bestsol.decs(i);
    cur_ground = bestsol.ground_list(i);
    ind_window = bestsol.index_window_guance(cur_satellite);
    
    t_s = bestsol.time_start_guance(cur_satellite);
    t_e = bestsol.time_end_guance(cur_satellite);
    
    t_s_window = Global.visible_window{cur_satellite,cur_ground}(2*bestsol.index_window_guance(cur_satellite)-1);
    t_e_window = Global.visible_window{cur_satellite,cur_ground}(2*bestsol.index_window_guance(cur_satellite));
    
    if t_s == 0 && t_e ==0
       continue; 
    end
    
    x = [t_s_window,t_e_window,t_e_window,t_s_window];
    y = [cur_ground-0.02,cur_ground-0.02,cur_ground+0.02,cur_ground+0.02]+i/50;
    patch(x,y,c_space(i,:),'facealpha',0.5);
    
    x = [t_s,t_e,t_e,t_s];
    y = [cur_ground-0.1,cur_ground-0.1,cur_ground+0.1,cur_ground+0.1]+i/50;

    patch(x,y,c_space(i,:),'facealpha',0.8);
    text(t_s+50,cur_ground+i/50,num2str(cur_satellite),'FontWeight','Bold','fontsize',16);
    ylim([0,5])
end
for i=0.5:1:4.5
   line([min(bestsol.time_start_guance),max(bestsol.time_end_guance)],[i,i],'linestyle','-.') 
end
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
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70