利用geatpy进化算法优化卫星观测覆盖面积

在卫星观测过程中，在指定的时间范围内，对于感兴趣区域的观测可能难以完整覆盖，但是部分区域又可能有冗余的观测，由于近年来的卫星多数都拥有摆动能力，所以如何利用卫星的摆动能力，使得对研究区的观测达到最优是一个要解决的问题。

由于网上现有的利用进化算法进行卫星覆盖面积优化的代码目前很少有开源的，这里只是提供一个优化的思路，观测条带的计算可以参考相关的论文。

 1、首先定义问题，这里就是单目标优化，计算最大覆盖面积

class MyProblem(ea.Problem): # 继承Problem父类
    def __init__(self):
        name = 'MyProblem' # 初始化name（函数名称，可以随意设置）
        M = 1 # 初始化M（目标维数）
        # 初始化maxormins（目标最小最大化标记列表，1：最小化；-1：最大化）
        maxormins = [-1]
        Dim = varnumber # 初始化Dim（决策变量维数,也就代表卫星条带的数量）
        # 初始化决策变量的类型，元素为0表示变量是连续的；1为离散的
        varTypes = [1] * varnumber # 初始化决策变量类型，0：连续；1：离散
        lb = [0] * varnumber# 决策变量下界
        ub = [1] * varnumber # 决策变量上界
        lbin = [1] * varnumber # 决策变量下边界
        ubin = [1] * varnumber# 决策变量上边界
        # 调用父类构造方法完成实例化
        ea.Problem.__init__(self, name, M, maxormins, Dim, varTypes, lb,ub, lbin, ubin)
    def aimFunc(self, pop): # 目标函数，pop为传入的种群对象
        
        Vars = pop.Phen # 得到决策变量矩阵
        areas=area_calculator3(Vars)# 计算目标函数值，赋值给pop种群对象的ObjV属性(面积)
        row_sums = areas.sum(axis=1)
        row_sums_matrix = row_sums.reshape(-1, 1)
        pop.ObjV = row_sums_matrix
        
        # 采用可行性法则处理约束，生成种群个体违反约束程度矩阵
        constraints=generate_constraints_based_on_ranges2(Vars, gdf)
        pop.CV = constraints

2、根据染色体和gdf数据计算染色体所代表的面积，gdf也就是存储卫星条带的geodataframe

def area_calculator3(selected): 
    nonzero_indices = np.argwhere(selected != 0)# 使用矢量化操作找出非零元素的索引
    computed_areas  = {}# 创建一个字典来存储已经计算过的行索引和相应的面积结果
    area_list = np.zeros(selected.shape[0])# 根据非零索引创建一个空的列表来存储面积
    
    # 遍历每一行的非零索引，并计算面积
    for row_index in np.unique(nonzero_indices[:, 0]):
        row_nonzero_indices = nonzero_indices[nonzero_indices[:, 0] == row_index, 1] # 获取当前行的非零元素索引
        
        # 将非零元素索引转换为不可变的元组，作为字典的键
        row_nonzero_indices_tuple = tuple(row_nonzero_indices)
        
        if row_nonzero_indices.size > 0:# 如果当前行有非零元素，则计算合并后的多边形与roi的交集面积

            if row_nonzero_indices_tuple in computed_areas:
                # 如果已计算过，直接从字典获取面积
                area_list[row_index] = computed_areas[row_nonzero_indices_tuple]
                
            else :
                polygons = [Polygon(gdf.iloc[i]['Rectangle']) for i in row_nonzero_indices]
                # 合并多边形
                merged_polygon = cascaded_union(polygons)
                # 计算多边形与 roi 的交集
                intersection = merged_polygon.intersection(roi)
                # 获取交集面积
                area = intersection.area
                area_list[row_index] = area
                computed_areas[row_nonzero_indices_tuple] = area
          
    # 将列表转换为二维数组以匹配原始函数的输出
    areaarray = area_list[:, np.newaxis]  
    return areaarray

3、定义约束条件，每个卫星的摆动角度在同一天中只能有1个

def generate_constraints_based_on_ranges2(Vars, gdf): #这个不适配多传感器
    # 基于日期、卫星和传感器进行分组
    grouped = gdf.groupby(['date', 'Satellite', 'sensor_name'])

    # 计算每组的索引范围
    index_ranges = grouped.apply(lambda x: (x.index.min(), x.index.max()))

    # 用于存储约束条件的列表
    constraints = []

    # 遍历每个卫星和传感器组合的索引范围
    for (_, _, _), (start_idx, end_idx) in index_ranges.items():
        if start_idx == end_idx:
            continue  # 如果只有一个观测机会，则跳过

        # 对每个卫星的每个传感器，生成摆动角度选择的约束条件
        columns = Vars[:, start_idx:end_idx+1]
        constraint = np.abs(np.sum(columns, axis=1) - 1)
        constraints.append(constraint)

    # 生成最终的约束条件数组
    constraints_array = np.vstack(constraints).T
    return constraints_array

4、问题优化

    # """============================实例化问题对象========================"""
        problem = MyProblem() # 实例化问题对象
        # """==============================种群设置==========================="""
        Encoding = 'RI' # 编码方式
        # Encoding = ['RI']*varnumber
        NIND = 40 # 种群规模
        Field = ea.crtfld(Encoding, problem.varTypes, problem.ranges,problem.borders) # 创建区域描述器
        # Fields = [Field1, Field2]
        # Fields=[Fields]*Satellitenumber
        
        population = ea.Population(Encoding, Field, NIND) #
        # print('实例化完成')
        # 实例化种群对象（此时种群还没被真正初始化，仅仅是生成一个种群对象）
        # """===========================算法参数设置=========================="""
        myAlgorithm = ea.soea_SEGA_templet(problem, population) #实例化一个算法模板对象
        myAlgorithm.MAXGEN = 100 # 最大遗传代数
        myAlgorithm.mutOper.Pm = 0.5 #  变异
        myAlgorithm.recOper.XOVR = 0.2 # 设置交叉概率
        myAlgorithm.drawing = 1 # 设置绘图方式
        # """=====================调用算法模板进行种群进化====================="""
        # [population1, population2] = myAlgorithm.run() # 执行算法模板
        # population.save() # 把最后一代种群的信息保存到文件中
        res = ea.optimize(myAlgorithm,
                              verbose=False,
                              drawing=1,
                              outputMsg=False,
                              drawLog=False,
                              saveFlag=False)
        binary_array =res['Vars']
        selected_rows = gdf[binary_array[0] == 1]   #这个就是最终的优化结果，要根据位置去gdf中找到对应的观测条带

由于我对于geatpy也是入门新手，这些代码只是能完成功能，可能存在效率不高的问题，希望大家在使用时能及时指出问题