用NSGA2算法解决一个两目标优化的VRP问题随机生成一百个客户点和一个仓库的坐标横纵坐标在0-50以内给定客户的配送需求和期望配送车辆到达时间给定配送车辆的工作时间载重约束和每辆车使用的固定成本目标函数是最小总成本和最大总的客户满意度成本包括车的固定成本和配送行驶成本配送行驶成本与距离成正比总的客户满意度为车辆实际到达时间和期望配送车辆到达时间的比较如果期望配送车辆到达时间小于车辆实际到达时间则

由于题目较为复杂，以下仅提供一个基本的VRP问题的NSGA2算法实现，并未包含全部题目要求。

首先，定义一个客户点的类：

class Customer:
    def __init__(self, id, x, y, demand, service_time, time_window):
        self.id = id
        self.x = x
        self.y = y
        self.demand = demand
        self.service_time = service_time
        self.time_window = time_window

其中，id为客户点的编号，x和y为客户点的坐标，demand为客户点的需求量，service_time为客户点的服务时间，time_window为客户点的时间窗口。

定义一个车辆的类：

class Vehicle:
    def __init__(self, capacity, fixed_cost, variable_cost, max_working_time):
        self.capacity = capacity
        self.fixed_cost = fixed_cost
        self.variable_cost = variable_cost
        self.max_working_time = max_working_time
        self.route = []
        self.load = 0
        self.total_distance = 0
        self.total_time = 0

其中，capacity为车辆的容量，fixed_cost为车辆的固定成本，variable_cost为车辆的配送成本，max_working_time为车辆的最大工作时间，route为车辆的路径，load为车辆的载货量，total_distance为车辆的总行驶距离，total_time为车辆的总工作时间。

定义一个VRP问题的类：

class VRP:
    def __init__(self, customers, depot, vehicles):
        self.customers = customers
        self.depot = depot
        self.vehicles = vehicles
    
    def evaluate(self, solution):
        # 计算一个解的目标函数值
        pass
    
    def crossover(self, solution1, solution2):
        # 交叉操作
        pass
    
    def mutation(self, solution):
        # 变异操作
        pass
    
    def select(self, population):
        # 选择操作
        pass
    
    def solve(self, max_generations, population_size):
        # 解决问题
        pass

其中，customers为一百个客户点的列表，depot为仓库的坐标，vehicles为配送车辆的列表，evaluate方法用于计算一个解的目标函数值，crossover方法用于进行交叉操作，mutation方法用于进行变异操作，select方法用于进行选择操作，solve方法用于解决问题。

首先，我们需要生成初始的种群。我们可以随机生成若干个解，每个解表示每辆车的路径，如下所示：

def generate_initial_population(self, population_size):
    population = []
    for i in range(population_size):
        solution = []
        remaining_customers = self.customers[:]
        for vehicle in self.vehicles:
            route = [self.depot]
            load = 0
            while remaining_customers and load + remaining_customers[0].demand <= vehicle.capacity:
                customer = remaining_customers.pop(0)
                route.append(customer)
                load += customer.demand
            route.append(self.depot)
            vehicle.route = route
            vehicle.load = load
        population.append(solution)
    return population

在生成初始种群后，我们需要计算每个解的目标函数值。这里我们以最小化总成本和最大化总客户满意度为目标函数。

def evaluate(self, solution):
    total_cost = 0
    total_satisfaction = 0
    for vehicle in self.vehicles:
        route = vehicle.route
        load = 0
        distance = 0
        time = 0
        for i in range(1, len(route)):
            distance += math.sqrt((route[i].x - route[i - 1].x) ** 2 + (route[i].y - route[i - 1].y) ** 2)
            time = max(time + distance, route[i].time_window[0]) + route[i].service_time
            if time > route[i].time_window[1]:
                satisfaction = -1
            else:
                satisfaction = 0.5
            total_satisfaction += satisfaction
        variable_cost = distance * vehicle.variable_cost
        total_cost += vehicle.fixed_cost + variable_cost
    return total_cost, total_satisfaction

其中，distance为车辆的行驶距离，time为车辆的到达时间，satisfaction为客户的满意度。

接下来，我们需要进行交叉和变异操作。这里我们采用交叉操作的方式是随机选择两个解，并将它们的路径进行交叉，如下所示：

def crossover(self, solution1, solution2):
    index1 = random.randint(0, len(solution1) - 1)
    index2 = random.randint(0, len(solution2) - 1)
    route1 = solution1[index1].route
    route2 = solution2[index2].route
    crossover_point1 = random.randint(1, len(route1) - 2)
    crossover_point2 = random.randint(1, len(route2) - 2)
    new_route1 = route1[:crossover_point1] + route2[crossover_point2:-1] + [self.depot]
    new_route2 = route2[:crossover_point2] + route1[crossover_point1:-1] + [self.depot]
    solution1[index1].route = new_route1
    solution2[index2].route = new_route2

变异操作的方式是随机选择一个解，并将它的路径中的两个客户点进行交换，如下所示：

def mutation(self, solution):
    index = random.randint(0, len(solution) - 1)
    route = solution[index].route
    mutation_point1 = random.randint(1, len(route) - 2)
    mutation_point2 = random.randint(1, len(route) - 2)
    route[mutation_point1], route[mutation_point2] = route[mutation_point2], route[mutation_point1]

最后，我们需要进行选择操作。这里采用NSGA2算法中的拥挤度算法来进行选择。

def select(self, population, population_size):
    fronts = self.fast_non_dominated_sort(population)
    new_population = []
    i = 0
    while len(new_population) + len(fronts[i]) <= population_size:
        crowding_distance = self.crowding_distance(fronts[i])
        for solution in fronts[i]:
            solution.crowding_distance = crowding_distance[solution]
            new_population.append(solution)
        i += 1
    if len(new_population) < population_size:
        crowding_distance = self.crowding_distance(fronts[i])
        fronts[i].sort(key=lambda solution: solution.crowding_distance, reverse=True)
        new_population += fronts[i][:population_size - len(new_population)]
    return new_population

其中，fast_non_dominated_sort方法用于将种群按照非支配排序，crowding_distance方法用于计算每个解的拥挤度。

最后，我们可以将以上方法整合到solve方法中，如下所示：

def solve(self, max_generations, population_size):
    population = self.generate_initial_population(population_size)
    for generation in range(max_generations):
        for solution in population:
            solution.cost, solution.satisfaction = self.evaluate(solution)
        fronts = self.fast_non_dominated_sort(population)
        for i in range(len(fronts)):
            self.crowding_distance(fronts[i])
        new_population = []
        i = 0
        while len(new_population) + len(fronts[i]) <= population_size:
            crowding_distance = self.crowding_distance(fronts[i])
            for solution in fronts[i]:
                solution.crowding_distance = crowding_distance[solution]
                new_population.append(solution)
            i += 1
        if len(new_population) < population_size:
            crowding_distance = self.crowding_distance(fronts[i])
            fronts[i].sort(key=lambda solution: solution.crowding_distance, reverse=True)
            new_population += fronts[i][:population_size - len(new_population)]
        population = new_population
        for i in range(population_size):
            if random.random() < 0.9:
                if random.random() < 0.5:
                    index1 = random.randint(0, population_size - 1)
                    index2 = random.randint(0, population_size - 1)
                    self.crossover(population[index1], population[index2])
                else:
                    index = random.randint(0, population_size - 1)
                    self.mutation(population[index])
    return fronts[0]

最后，我们可以调用solve方法，求解VRP问题的帕累托前沿：

customers = []
for i in range(100):
    x = random.uniform(0, 50)
    y = random.uniform(0, 50)
    demand = random.randint(1, 10)
    service_time = demand
    time_window = (random.uniform(0, 10), random.uniform(20, 30))
    customers.append(Customer(i + 1, x, y, demand, service_time, time_window))
depot = (25, 25)
vehicles = [Vehicle(50, 1000, 2, 8) for i in range(10)]
vrp = VRP(customers, depot, vehicles)
fronts = vrp.solve(100, 100)
for solution in fronts[0]:
    print(solution.cost, solution.satisfaction)