运行下面的代码生成了一个仓库物资数量不动一直是0和征地的物资数量一直以每天的无补给的情况进行消耗减少import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import random
import math
from collections import defaultdict
import time
# 配置参数
POPULATION_SIZE = 50
GENERATIONS = 100
MUTATION_RATE = 0.15
ELITE_RATE = 0.1
CROSSOVER_RATE = 0.85
# 初始化数据
file_path = "C:\\Users\\Dd\\Desktop\\工作簿1.xlsx"
sheet_name = 'Sheet3'
df = pd.read_excel(file_path, sheet_name=sheet_name)
a = [tuple(row) for row in df[["栅格x坐标", "栅格y坐标", "高程(米)"]].values]
#z:[(10996, 6648, 500), (9979, 5931, 300), (9564, 4823, 450), (8570, 4707, 280), (8154, 3875, 320), (7739, 3390, 400), (6907, 2189, 210), (6560, 1173, 330), (7054, 410, 190)]
#c:[(5333, 8223, 280), (8039, 7941, 225), (5075, 6987, 180), (4698, 6162, 290), (3279, 5169, 120), (5922, 4615, 230), (4305, 3413, 175), (3265, 2674, 300), (4166, 760, 240)]
#h:[(3150, 10713, 100)]
location_z = a[0:9]
name_z = ['Z1', 'Z2', 'Z3', 'Z4', 'Z5', 'Z6', 'Z7', 'Z8', 'Z9']
location_c = a[9:18]
name_c = ['C1', 'C2', 'C3', 'C4', 'C5', 'C6', 'C7', 'C8', 'C9']
location_h = [a[18]]
name_h = ['H']
# 参数设置
one_day = 7.5
z_people = [68, 60, 38, 60, 90, 60, 38, 90, 60]
z_min = [people * one_day * 0.25 for people in z_people] # 最低库存(1/4日需求量)
z_max = [people * one_day * 3.5 for people in z_people] # 最大库存(3.5日需求量)
z_start = [1020.0, 900.0, 570.0, 900.0, 1350.0, 900.0, 570.0, 1350.0, 900.0]
a_time = 0.5
b_time = 1
c_max = [600, 800, 800, 800, 800, 800, 600, 800, 500]
a_m = 400
b_m = 600
a_v = 60 # km/h
b_v = 50 # km/h
a_consumption = 0.01 # 每公里耗电
b_consumption = 0.01 # 每公里耗电
DAYS = 15
WORK_HOURS = 24
# 改进的距离计算(考虑高程影响)
def calculate_distance(loc1, loc2):
dx = (loc1[0] - loc2[0]) * 5 # 转换为米
dy = (loc1[1] - loc2[1]) * 5
dz = loc1[2] - loc2[2]
return math.sqrt(dx ** 2 + dy ** 2 + dz ** 2) / 1000 # 转换为km
# 初始化距离矩阵(预计算)
distances = {
'c_to_h': [calculate_distance(loc, location_h[0]) for loc in location_c],
'c_to_z': [[calculate_distance(c, z) for z in location_z] for c in location_c],
'z_to_z': [[calculate_distance(z1, z2) for z2 in location_z] for z1 in location_z],
'z_to_c': [[calculate_distance(z, c) for c in location_c] for z in location_z]
}
# 辅助函数:充电
def charge_battery(current, hours):
"""充电函数,每小时恢复20%电量"""
return min(1.0, current + 0.2 * hours)
# 改进的个体表示
def create_individual():
ck_active = [random.random() < 0.8 for _ in range(9)]
if not any(ck_active):
ck_active[random.randint(0, 8)] = True
# 更合理的飞机数量范围
b_count = random.randint(8, 15)
b_locations = random.choices(name_c, k=b_count)
a_count = random.randint(8, 15)
a_locations = random.choices(name_c, k=a_count)
return {
'warehouses': ck_active,
'b_planes': {i: {'loc': loc, 'battery': 1.0, 'busy_until': 0}
for i, loc in enumerate(b_locations)},
'a_planes': {i: {'loc': loc, 'battery': 1.0, 'busy_until': 0}
for i, loc in enumerate(a_locations)}
}
def fitness(individual):
ck_active = individual['warehouses']
b_planes = individual['b_planes']
a_planes = individual['a_planes']
z_inventory = z_start.copy()
c_inventory = [0] * 9
penalty = 0
total_delivered = 0 # 总配送量
unmet_demand = 0 # 未满足需求
# 记录配送日志(仅用于调试)
delivery_log = []
# 简化调度逻辑:固定时间步长
for day in range(DAYS):
# 每日消耗
daily_consumption = [p * one_day for p in z_people]
for i in range(9):
z_inventory[i] -= daily_consumption[i]
if z_inventory[i] < z_min[i]:
penalty += (z_min[i] - z_inventory[i]) * 100
unmet_demand += z_min[i] - z_inventory[i]
if z_inventory[i] < 0:
penalty += abs(z_inventory[i]) * 1000 # 严重短缺惩罚
# 模拟每小时
for hour in range(WORK_HOURS):
current_time = day * WORK_HOURS + hour
# B飞机调度 (从H到C) - 简化逻辑
for pid, plane in list(b_planes.items()):
if current_time >= plane['busy_until']:
# 飞机在H点或C点
if plane['loc'] == 'H':
# 寻找库存空间最大的仓库
best_c = None
max_capacity = -1
for ci, active in enumerate(ck_active):
if not active: continue
capacity = c_max[ci] - c_inventory[ci]
if capacity > max_capacity:
max_capacity = capacity
best_c = ci
if best_c is not None and max_capacity > 0:
dist = distances['c_to_h'][best_c]
flight_time = dist / b_v
consumption = dist * b_consumption
if plane['battery'] >= consumption:
# 计算可运输量
cargo = min(b_m, max_capacity)
# 执行运输
c_inventory[best_c] += cargo
total_delivered += cargo
# 更新飞机状态
plane['busy_until'] = current_time + flight_time + b_time
plane['battery'] -= consumption
plane['loc'] = name_c[best_c]
# 记录配送
delivery_log.append({
'time': current_time,
'type': 'B',
'from': 'H',
'to': name_c[best_c],
'cargo': cargo,
'battery': plane['battery']
})
else:
# 电量不足,充电
plane['busy_until'] = current_time + 1
plane['battery'] = charge_battery(plane['battery'], 1)
else:
# 没有可用仓库,原地待命
plane['busy_until'] = current_time + 1
else:
# 飞机在C点,返回H点
c_idx = name_c.index(plane['loc'])
dist = distances['c_to_h'][c_idx]
flight_time = dist / b_v
consumption = dist * b_consumption
if plane['battery'] >= consumption:
plane['busy_until'] = current_time + flight_time + b_time
plane['battery'] -= consumption
plane['loc'] = 'H'
else:
# 电量不足,充电
plane['busy_until'] = current_time + 1
plane['battery'] = charge_battery(plane['battery'], 1)
# A飞机调度 (从C到Z) - 优化版本
for pid, plane in list(a_planes.items()):
if current_time >= plane['busy_until']:
# 飞机在C点 - 准备向Z点运输
if plane['loc'] in name_c:
c_idx = name_c.index(plane['loc'])
# 检查仓库是否启用且有足够库存
if not ck_active[c_idx] or c_inventory[c_idx] < a_m:
# 条件不满足,充电或待命
plane['busy_until'] = current_time + 1
plane['battery'] = min(1.0, plane['battery'] + 0.2)
continue
# 寻找最需要补给的Z点(综合考虑库存和距离)
best_z = None
best_score = -float('inf')
for zi in range(9):
# 计算Z点需求紧急程度 (0-1)
urgency = 1 - min(1.0, z_inventory[zi] / z_min[zi]) if z_min[zi] > 0 else 1
# 计算距离因素 (0-1)
distance = distances['c_to_z'][c_idx][zi]
distance_factor = 1 - min(1.0, distance / 50) # 假设最大距离50km
# 综合评分 (紧急程度70%,距离30%)
score = 0.7 * urgency + 0.3 * distance_factor
# 只考虑有需求且未满的Z点
if z_inventory[zi] < z_max[zi] and score > best_score:
best_score = score
best_z = zi
if best_z is not None:
dist = distances['c_to_z'][c_idx][best_z]
flight_time_needed = dist / a_v
consumption = dist * a_consumption
# 计算可运输量 (考虑仓库库存和Z点容量)
cargo = min(
a_m,
c_inventory[c_idx],
z_max[best_z] - z_inventory[best_z]
)
# 检查电量是否足够单程飞行
if plane['battery'] >= consumption and cargo > 0:
# 执行运输
z_inventory[best_z] += cargo # 更新Z点库存
c_inventory[c_idx] -= cargo # 更新仓库库存
total_delivered += cargo # 更新总配送量
# 更新飞机状态
plane['busy_until'] = current_time + flight_time_needed + a_time
plane['battery'] -= consumption
plane['loc'] = name_z[best_z] # 位置更新到Z点
# 记录配送日志
delivery_log.append({
'time': current_time,
'type': 'A',
'from': name_c[c_idx],
'to': name_z[best_z],
'cargo': cargo,
'battery': plane['battery']
})
else:
# 电量不足或无法运输,充电
plane['busy_until'] = current_time + 1
plane['battery'] = min(1.0, plane['battery'] + 0.2)
else:
# 没有需要补给的Z点
plane['busy_until'] = current_time + 1
# 飞机在Z点 - 需要返回C点
elif plane['loc'] in name_z:
z_idx = name_z.index(plane['loc'])
# 寻找最近的启用仓库
best_c = None
min_distance = float('inf')
for ci in range(9):
if ck_active[ci]:
dist = distances['z_to_c'][z_idx][ci]
if dist < min_distance:
min_distance = dist
best_c = ci
if best_c is not None:
dist = min_distance
flight_time_needed = dist / a_v
consumption = dist * a_consumption
# 检查电量是否足够返回
if plane['battery'] >= consumption:
# 执行返回
plane['busy_until'] = current_time + flight_time_needed + a_time
plane['battery'] -= consumption
plane['loc'] = name_c[best_c] # 位置更新到C点
# 记录返回日志
delivery_log.append({
'time': current_time,
'type': 'A_RETURN',
'from': name_z[z_idx],
'to': name_c[best_c],
'cargo': 0,
'battery': plane['battery']
})
else:
# 电量不足,充电
plane['busy_until'] = current_time + 1
plane['battery'] = min(1.0, plane['battery'] + 0.2)
else:
# 没有可用仓库,原地待命
plane['busy_until'] = current_time + 1
# 每小时检查仓库容量
for ci in range(9):
if ck_active[ci] and c_inventory[ci] > c_max[ci]:
penalty += (c_inventory[ci] - c_max[ci]) * 10
c_inventory[ci] = c_max[ci]
# 最终库存检查
for i in range(9):
if z_inventory[i] < z_min[i]:
penalty += (z_min[i] - z_inventory[i]) * 100
# 计算适应度:奖励配送量,惩罚短缺和未满足需求
fitness_value = total_delivered - penalty - unmet_demand
# 确保适应度非负
return max(0.001, fitness_value)
# 改进的选择、交叉和变异
def tournament_selection(population, fitnesses, k=3):
selected = []
for _ in range(len(population)):
candidates = random.sample(list(zip(population, fitnesses)), k)
winner = max(candidates, key=lambda x: x[1])[0]
selected.append(winner)
return selected
def crossover(parent1, parent2):
if random.random() < CROSSOVER_RATE:
# 仓库状态交叉
crossover_point = random.randint(1, 8)
new_warehouses = parent1['warehouses'][:crossover_point] + parent2['warehouses'][crossover_point:]
# 飞机交叉:随机选择父代的飞机
new_b_planes = {}
b_planes1 = list(parent1['b_planes'].values())
b_planes2 = list(parent2['b_planes'].values())
b_count = random.randint(min(len(b_planes1), len(b_planes2)), max(len(b_planes1), len(b_planes2)))
for i in range(b_count):
if i < len(b_planes1) and random.random() < 0.5:
new_b_planes[i] = b_planes1[i].copy()
elif i < len(b_planes2):
new_b_planes[i] = b_planes2[i].copy()
new_a_planes = {}
a_planes1 = list(parent1['a_planes'].values())
a_planes2 = list(parent2['a_planes'].values())
a_count = random.randint(min(len(a_planes1), len(a_planes2)), max(len(a_planes1), len(a_planes2)))
for i in range(a_count):
if i < len(a_planes1) and random.random() < 0.5:
new_a_planes[i] = a_planes1[i].copy()
elif i < len(a_planes2):
new_a_planes[i] = a_planes2[i].copy()
return {
'warehouses': new_warehouses,
'b_planes': new_b_planes,
'a_planes': new_a_planes
}
return parent1.copy()
def mutate(individual):
# 仓库变异
for i in range(9):
if random.random() < MUTATION_RATE:
individual['warehouses'][i] = not individual['warehouses'][i]
if not any(individual['warehouses']):
individual['warehouses'][random.randint(0, 8)] = True
# B飞机变异
for pid, plane in list(individual['b_planes'].items()):
if random.random() < MUTATION_RATE:
if random.random() < 0.3 and len(individual['b_planes']) > 5:
del individual['b_planes'][pid]
else:
plane['loc'] = random.choice(name_c)
# 随机添加B飞机
if random.random() < MUTATION_RATE * 0.5 and len(individual['b_planes']) < 15:
new_pid = max(individual['b_planes'].keys()) + 1 if individual['b_planes'] else 0
individual['b_planes'][new_pid] = {
'loc': random.choice(name_c),
'battery': min(1.0, random.uniform(0.5, 1.0)),
'busy_until': 0
}
# A飞机变异
for pid, plane in list(individual['a_planes'].items()):
if random.random() < MUTATION_RATE:
if random.random() < 0.3 and len(individual['a_planes']) > 5:
del individual['a_planes'][pid]
else:
plane['loc'] = random.choice(name_c)
# 随机添加A飞机
if random.random() < MUTATION_RATE * 0.5 and len(individual['a_planes']) < 15:
new_pid = max(individual['a_planes'].keys()) + 1 if individual['a_planes'] else 0
individual['a_planes'][new_pid] = {
'loc': random.choice(name_c),
'battery': min(1.0, random.uniform(0.5, 1.0)),
'busy_until': 0
}
return individual
# 改进的遗传算法
def genetic_algorithm():
population = [create_individual() for _ in range(POPULATION_SIZE)]
best_fitness = -float('inf')
best_individual = None
fitness_history = []
start_time = time.time()
for gen in range(GENERATIONS):
# 评估适应度
fitnesses = []
for ind in population:
fit = fitness(ind)
fitnesses.append(fit)
# 更新最佳解
current_best = max(fitnesses)
current_idx = fitnesses.index(current_best)
if current_best > best_fitness:
best_fitness = current_best
best_individual = population[current_idx].copy()
print(f"Generation {gen + 1}: New best fitness = {best_fitness:.2f}")
fitness_history.append(best_fitness)
# 选择
selected = tournament_selection(population, fitnesses)
# 交叉
children = []
for i in range(0, len(selected), 2):
if i + 1 < len(selected):
child1 = crossover(selected[i], selected[i + 1])
child2 = crossover(selected[i + 1], selected[i])
children.extend([child1, child2])
else:
children.append(selected[i].copy())
# 变异
next_population = [mutate(ind) for ind in children]
# 精英保留
elite_size = int(POPULATION_SIZE * ELITE_RATE)
if elite_size > 0:
elite_indices = np.argsort(fitnesses)[-elite_size:]
for idx in elite_indices:
next_population.append(population[idx].copy())
# 保持种群大小
population = next_population[:POPULATION_SIZE]
end_time = time.time()
print(f"\nGenetic algorithm completed in {end_time - start_time:.2f} seconds")
return best_individual, best_fitness, fitness_history
# 运行算法
best_solution, best_fitness, history = genetic_algorithm()
# 结果分析
print("\n=== Best Solution ===")
print(f"Fitness: {best_fitness:.2f}")
print(f"Active Warehouses: {[name_c[i] for i, active in enumerate(best_solution['warehouses']) if active]}")
print(f"B Planes: {len(best_solution['b_planes'])}")
print(f"A Planes: {len(best_solution['a_planes'])}")
# 可视化
plt.figure(figsize=(12, 6))
plt.plot(history, 'b-', linewidth=2)
plt.title('Fitness Progression')
plt.xlabel('Generation')
plt.ylabel('Fitness')
plt.grid(True, linestyle='--', alpha=0.7)
plt.tight_layout()
plt.savefig('fitness_progression.png', dpi=300)
plt.show()
# 模拟最佳解决方案的物流情况
def simulate_best_solution(solution):
print("\nSimulating best solution...")
ck_active = solution['warehouses']
b_planes = solution['b_planes']
a_planes = solution['a_planes']
z_inventory = z_start.copy()
c_inventory = [0] * 9
total_delivered = 0
# 每日库存记录
daily_z_inventory = []
daily_c_inventory = []
for day in range(DAYS):
# 每日消耗
daily_consumption = [p * one_day for p in z_people]
for i in range(9):
z_inventory[i] -= daily_consumption[i]
# 记录每日库存
daily_z_inventory.append(z_inventory.copy())
daily_c_inventory.append(c_inventory.copy())
# 模拟每小时
for hour in range(WORK_HOURS):
current_time = day * WORK_HOURS + hour
# B飞机调度
for pid, plane in b_planes.items():
if current_time >= plane['busy_until']:
if plane['loc'] == 'H':
for ci, active in enumerate(ck_active):
if active and c_inventory[ci] < c_max[ci]:
dist = distances['c_to_h'][ci]
flight_time = dist / b_v
consumption = dist * b_consumption
if plane['battery'] >= consumption:
cargo = min(b_m, c_max[ci] - c_inventory[ci])
c_inventory[ci] += cargo
total_delivered += cargo
plane['busy_until'] = current_time + flight_time + b_time
plane['battery'] -= consumption
plane['loc'] = name_c[ci]
print(f"Day {day + 1} Hour {hour}: B{pid} delivered {cargo}kg to {name_c[ci]}")
else:
c_idx = name_c.index(plane['loc'])
dist = distances['c_to_h'][c_idx]
flight_time = dist / b_v
consumption = dist * b_consumption
if plane['battery'] >= consumption:
plane['busy_until'] = current_time + flight_time + b_time
plane['battery'] -= consumption
plane['loc'] = 'H'
# A飞机调度
for pid, plane in a_planes.items():
if current_time >= plane['busy_until']:
if plane['loc'] in name_c:
c_idx = name_c.index(plane['loc'])
if ck_active[c_idx] and c_inventory[c_idx] > 0:
# 寻找库存最低的Z点
best_z = None
min_inventory = float('inf')
for zi in range(9):
if z_inventory[zi] < min_inventory and z_inventory[zi] < z_max[zi]:
min_inventory = z_inventory[zi]
best_z = zi
if best_z is not None:
dist = distances['c_to_z'][c_idx][best_z]
flight_time = dist / a_v
consumption = dist * a_consumption
if plane['battery'] >= consumption:
cargo = min(a_m, c_inventory[c_idx], z_max[best_z] - z_inventory[best_z])
z_inventory[best_z] += cargo
c_inventory[c_idx] -= cargo
total_delivered += cargo
plane['busy_until'] = current_time + flight_time + a_time
plane['battery'] -= consumption
plane['loc'] = name_z[best_z]
print(f"Day {day + 1} Hour {hour}: A{pid} delivered {cargo}kg to {name_z[best_z]}")
print(f"\nTotal delivered: {total_delivered}kg")
return daily_z_inventory, daily_c_inventory
# 模拟最佳解决方案
daily_z_inv, daily_c_inv = simulate_best_solution(best_solution)
# 可视化库存变化
plt.figure(figsize=(14, 10))
plt.suptitle('Inventory Levels Over Time', fontsize=16)
# Z点库存
plt.subplot(2, 1, 1)
for i in range(9):
inventory = [day[i] for day in daily_z_inv]
plt.plot(range(1, DAYS + 1), inventory, label=f'{name_z[i]}', marker='o')
plt.axhline(y=0, color='r', linestyle='--', alpha=0.3)
plt.xlabel('Day')
plt.ylabel('Inventory (kg)')
plt.title('Frontline Base (Z) Inventory')
plt.legend()
plt.grid(True)
# C点库存
plt.subplot(2, 1, 2)
for i in range(9):
if best_solution['warehouses'][i]:
inventory = [day[i] for day in daily_c_inv]
plt.plot(range(1, DAYS + 1), inventory, label=f'{name_c[i]}', marker='s')
plt.xlabel('Day')
plt.ylabel('Inventory (kg)')
plt.title('Warehouse (C) Inventory')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.savefig('inventory_levels.png', dpi=300)
plt.show()
最新发布