# 02.4 空间优化 ## 核心问题 > 如何在无穷可能中找到"最优"的空间配置？ > 当目标相互冲突时，什么是可接受的妥协解？ --- ## 概念讲解 ### 什么是空间优化 空间优化是在空间约束下寻找最优决策方案的过程： ``` ┌─────────────────────────────────────────────────────────────┐ │ 空间优化问题的结构 │ ├─────────────────────────────────────────────────────────────┤ │ │ │ ┌─────────────────────────────────────────────────────┐ │ │ │ 目标函数 │ │ │ │ Objective Function = f(x, y, ...) │ │ │ │ │ │ │ │ 例: 最大化生态连通性 │ │ │ │ Maximize Σ connectivity(patch_i, patch_j) │ │ │ │ │ │ │ │ 可能是: │ │ │ │ - 单目标优化 (一个目标) │ │ │ │ - 多目标优化 (多个目标，需权衡) │ │ │ └─────────────────────────────────────────────────────┘ │ │ │ │ │ ▼ │ │ ┌─────────────────────────────────────────────────────┐ │ │ │ 决策变量 │ │ │ │ Decision Variables = X │ │ │ │ │ │ │ │ 例: 哪些位置建立生态廊道 │ │ │ │ X = [0, 1, 0, 1, 1, ...] │ │ │ │ (1=建设, 0=不建设) │ │ │ └─────────────────────────────────────────────────────┘ │ │ │ │ │ ▼ │ │ ┌─────────────────────────────────────────────────────┐ │ │ │ 约束条件 │ │ │ │ Constraints = g(X) ≤ 0 │ │ │ │ │ │ │ │ 例: │ │ │ │ - 预算约束: Σ cost ≤ budget │ │ │ │ - 空间约束: 避开建设用地 │ │ │ │ - 连通性约束: 每个源地至少连接一条廊道 │ │ │ └─────────────────────────────────────────────────────┘ │ │ │ │ │ ▼ │ │ ┌─────────────────────────────────────────────────────┐ │ │ │ 可行域 │ │ │ │ Feasible Region │ │ │ │ │ │ │ │ 满足所有约束的解空间 │ │ │ │ 在可行域内寻找使目标函数最优的解 │ │ │ └─────────────────────────────────────────────────────┘ │ │ │ └─────────────────────────────────────────────────────────────┘ ``` ### 空间优化问题分类 | 问题类型 | 目标 | 典型应用 | 求解难度 | |---------|------|---------|---------| | **选址问题** | 选定最优位置 | 设施选址、生态源地识别 | NP-hard | | **覆盖问题** | 覆盖最大需求 | 保护区设计、服务覆盖 | NP-hard | | **分配问题** | 最优分配资源 | 土地利用分配 | NP-hard | | **路径问题** | 最短/最优路径 | 廊道设计、路线规划 | P (单点对点) | | **布局问题** | 优化空间布局 | 城市规划、景观设计 | NP-hard | | **网络设计** | 优化网络结构 | 生态网络、交通网络 | NP-hard | ### 求解方法谱系 ``` ┌─────────────────────────────────────────────────────────────┐ │ 优化求解方法 │ ├─────────────────────────────────────────────────────────────┤ │ │ │ 1. 精确算法 (Exact Methods) │ │ ┌─────────────────────────────────────────────────┐ │ │ │ - 线性规划 (LP): 单纯形法、内点法 │ │ │ │ - 整数规划 (IP): 分支定界、割平面 │ │ │ │ - 动态规划 (DP): 最优子结构 │ │ │ │ │ │ │ │ 优点: 保证全局最优 │ │ │ │ 缺点: 只适用于小规模问题 │ │ │ └─────────────────────────────────────────────────┘ │ │ │ │ 2. 启发式算法 (Heuristics) │ │ ┌─────────────────────────────────────────────────┐ │ │ │ - 贪心算法: 每步选择局部最优 │ │ │ │ - 构造式算法: 逐步构建解 │ │ │ │ │ │ │ │ 优点: 快速、简单 │ │ │ │ 缺点: 不保证最优 │ │ │ └─────────────────────────────────────────────────┘ │ │ │ │ 3. 元启发式算法 (Metaheuristics) │ │ ┌─────────────────────────────────────────────────┐ │ │ │ - 遗传算法 (GA): 模拟进化 │ │ │ │ - 模拟退火 (SA): 模拟金属退火 │ │ │ │ - 蚁群算法 (ACO): 模拟蚂蚁觅食 │ │ │ │ - 粒子群优化 (PSO): 模拟鸟群 │ │ │ │ │ │ │ │ 优点: 可处理大规模、非线性问题 │ │ │ │ 缺点: 参数敏感，不保证全局最优 │ │ │ └─────────────────────────────────────────────────┘ │ │ │ └─────────────────────────────────────────────────────────────┘ ``` ### 常见空间优化模型 #### 1. p-中值问题 (p-Median Problem) 选择p个设施，使所有需求点到最近设施的距离之和最小。 ``` Minimize: Σ Σ demand_i × distance(i, j) × x(i,j) i j∈selected Subject to: - 选恰好p个设施: Σ y_j = p j - 每个需求点被服务: Σ x(i,j) = 1, ∀i j - 只有被选中的设施才能服务: x(i,j) ≤ y_j, ∀i,j ``` #### 2. 最大覆盖问题 (Maximal Covering Problem) 用p个设施覆盖尽可能多的需求。 ``` Maximize: Σ demand_i × y_i i Subject to: - 选恰好p个设施: Σ x_j = p j - 覆盖关系: y_i ≤ Σ x_j, ∀i (j在i的覆盖范围内) j∈N(i) - 选恰好p个: Σ x_j = p ``` #### 3. 生态廊道优化 在预算约束下最大化生态连通性。 ``` Maximize: Σ connectivity_gain(c) × x_c c∈candidates Subject to: - 预算约束: Σ cost(c) × x_c ≤ budget c∈candidates - 连通性约束: 每个源地至少有一条廊道连接 ``` --- ## 设计原理 ### 遗传算法 (Genetic Algorithm) 遗传算法模拟自然进化过程，是空间优化中最常用的元启发式方法： ```python import numpy as np from typing import Callable, List, Tuple, Optional, Dict import random class GeneticAlgorithm: """ 遗传算法实现 核心思想：模拟自然选择、交叉、变异 """ def __init__(self, objective_func: Callable, n_variables: int, variable_type: str = 'binary', bounds: Optional[Tuple] = None, population_size: int = 100, mutation_rate: float = 0.01, crossover_rate: float = 0.8, elite_size: int = 2): """ Args: objective_func: 目标函数 (最小化) n_variables: 决策变量数量 variable_type: 'binary' 或 'continuous' bounds: (min, max) 连续变量的边界 population_size: 种群大小 mutation_rate: 变异率 crossover_rate: 交叉率 elite_size: 精英保留数量 """ self.objective_func = objective_func self.n_variables = n_variables self.variable_type = variable_type self.bounds = bounds or (0, 1) self.population_size = population_size self.mutation_rate = mutation_rate self.crossover_rate = crossover_rate self.elite_size = elite_size self.population = None self.fitness = None self.best_solution = None self.best_fitness = float('inf') self.history = [] def initialize(self): """初始化种群""" if self.variable_type == 'binary': self.population = np.random.randint( 0, 2, (self.population_size, self.n_variables) ) else: # continuous self.population = np.random.uniform( self.bounds[0], self.bounds[1], (self.population_size, self.n_variables) ) def evaluate(self): """评估种群适应度""" self.fitness = np.array([ self.objective_func(individual) for individual in self.population ]) # 更新最优解 best_idx = np.argmin(self.fitness) if self.fitness[best_idx] < self.best_fitness: self.best_fitness = self.fitness[best_idx] self.best_solution = self.population[best_idx].copy() self.history.append(self.best_fitness) def selection(self, method: str = 'tournament') -> np.ndarray: """ 选择操作 Args: method: 'tournament'(锦标赛) 或 'roulette'(轮盘赌) """ selected = [] if method == 'tournament': tournament_size = 3 for _ in range(self.population_size - self.elite_size): # 随机选择tournament_size个个体 candidates = np.random.choice( self.population_size, tournament_size, replace=False ) # 选择适应度最好的 winner = candidates[np.argmin(self.fitness[candidates])] selected.append(self.population[winner].copy()) elif method == 'roulette': # 转换为适应度（越小越好→越大越好） fitness_values = self.fitness if fitness_values.min() < 0: fitness_values = fitness_values - fitness_values.min() + 1 # 归一化 probs = 1 / fitness_values probs = probs / probs.sum() for _ in range(self.population_size - self.elite_size): idx = np.random.choice(self.population_size, p=probs) selected.append(self.population[idx].copy()) return np.array(selected) def crossover(self, parent1: np.ndarray, parent2: np.ndarray) -> Tuple: """ 交叉操作 Args: parent1, parent2: 父代个体 Returns: 两个子代个体 """ if self.variable_type == 'binary': # 单点交叉 if np.random.random() < self.crossover_rate: point = np.random.randint(1, self.n_variables) child1 = np.concatenate([parent1[:point], parent2[point:]]) child2 = np.concatenate([parent2[:point], parent1[point:]]) else: child1, child2 = parent1.copy(), parent2.copy() else: # continuous # 模拟二进制交叉 (SBX) if np.random.random() < self.crossover_rate: eta = 2 # 分布指数 u = np.random.random(self.n_variables) beta = np.where( u <= 0.5, (2 * u) ** (1 / (eta + 1)), (1 / (2 * (1 - u))) ** (1 / (eta + 1)) ) child1 = 0.5 * ((1 + beta) * parent1 + (1 - beta) * parent2) child2 = 0.5 * ((1 - beta) * parent1 + (1 + beta) * parent2) # 边界处理 child1 = np.clip(child1, self.bounds[0], self.bounds[1]) child2 = np.clip(child2, self.bounds[0], self.bounds[1]) else: child1, child2 = parent1.copy(), parent2.copy() return child1, child2 def mutate(self, individual: np.ndarray) -> np.ndarray: """ 变异操作 Args: individual: 个体 Returns: 变异后的个体 """ mutated = individual.copy() if self.variable_type == 'binary': # 位翻转变异 mask = np.random.random(self.n_variables) < self.mutation_rate mutated[mask] = 1 - mutated[mask] else: # continuous # 多项式变异 for i in range(self.n_variables): if np.random.random() < self.mutation_rate: delta = np.random.normal(0, 0.1 * (self.bounds[1] - self.bounds[0])) mutated[i] = np.clip( mutated[i] + delta, self.bounds[0], self.bounds[1] ) return mutated def evolve(self, n_generations: int) -> Dict: """ 进化指定代数 Args: n_generations: 进化代数 Returns: 结果字典 """ self.initialize() for generation in range(n_generations): # 评估 self.evaluate() # 精英保留 elite_indices = np.argsort(self.fitness)[:self.elite_size] elite = self.population[elite_indices].copy() # 选择 selected = self.selection(method='tournament') # 交叉 offspring = [] for i in range(0, len(selected), 2): if i + 1 < len(selected): child1, child2 = self.crossover(selected[i], selected[i+1]) offspring.extend([child1, child2]) else: offspring.append(selected[i]) # 变异 offspring = np.array([self.mutate(ind) for ind in offspring]) # 组合精英和后代 self.population = np.vstack([elite, offspring]) # 确保种群大小 if len(self.population) > self.population_size: self.population = self.population[:self.population_size] return { 'best_solution': self.best_solution, 'best_fitness': self.best_fitness, 'history': self.history } ``` ### 模拟退火算法 (Simulated Annealing) 模拟退火模拟金属冷却过程，能跳出局部最优： ```python import numpy as np from typing import Callable, Tuple, Optional import math class SimulatedAnnealing: """ 模拟退火算法 核心思想：以概率接受劣解，避免陷入局部最优 """ def __init__(self, objective_func: Callable, n_variables: int, variable_type: str = 'binary', bounds: Optional[Tuple] = None, initial_temp: float = 1000, cooling_rate: float = 0.95, min_temp: float = 0.01): """ Args: objective_func: 目标函数 (最小化) n_variables: 决策变量数量 variable_type: 'binary' 或 'continuous' bounds: 连续变量边界 initial_temp: 初始温度 cooling_rate: 降温率 min_temp: 最低温度 """ self.objective_func = objective_func self.n_variables = n_variables self.variable_type = variable_type self.bounds = bounds or (0, 1) self.initial_temp = initial_temp self.cooling_rate = cooling_rate self.min_temp = min_temp self.current_solution = None self.current_fitness = None self.best_solution = None self.best_fitness = float('inf') self.history = [] def initialize(self) -> np.ndarray: """生成初始解""" if self.variable_type == 'binary': return np.random.randint(0, 2, self.n_variables) else: return np.random.uniform( self.bounds[0], self.bounds[1], self.n_variables ) def generate_neighbor(self, solution: np.ndarray) -> np.ndarray: """生成邻域解""" neighbor = solution.copy() if self.variable_type == 'binary': # 随机翻转一位 idx = np.random.randint(self.n_variables) neighbor[idx] = 1 - neighbor[idx] else: # continuous # 在随机维度添加小扰动 idx = np.random.randint(self.n_variables) delta = np.random.normal(0, 0.1 * (self.bounds[1] - self.bounds[0])) neighbor[idx] = np.clip( neighbor[idx] + delta, self.bounds[0], self.bounds[1] ) return neighbor def accept_probability(self, current_fitness: float, new_fitness: float, temperature: float) -> float: """ 计算接受概率 (Metropolis准则) Args: current_fitness: 当前解适应度 new_fitness: 新解适应度 temperature: 当前温度 Returns: 接受概率 """ if new_fitness < current_fitness: return 1.0 # 更优解，一定接受 else: # 劣解以概率接受 return math.exp(-(new_fitness - current_fitness) / temperature) def optimize(self, max_iterations: int = 10000) -> Dict: """ 执行优化 Args: max_iterations: 最大迭代次数 Returns: 结果字典 """ # 初始化 self.current_solution = self.initialize() self.current_fitness = self.objective_func(self.current_solution) self.best_solution = self.current_solution.copy() self.best_fitness = self.current_fitness temperature = self.initial_temp for iteration in range(max_iterations): # 生成邻域解 neighbor = self.generate_neighbor(self.current_solution) neighbor_fitness = self.objective_func(neighbor) # 决定是否接受 prob = self.accept_probability( self.current_fitness, neighbor_fitness, temperature ) if np.random.random() < prob: self.current_solution = neighbor self.current_fitness = neighbor_fitness # 更新最优解 if self.current_fitness < self.best_fitness: self.best_fitness = self.current_fitness self.best_solution = self.current_solution.copy() # 记录 self.history.append(self.best_fitness) # 降温 temperature *= self.cooling_rate if temperature < self.min_temp: break return { 'best_solution': self.best_solution, 'best_fitness': self.best_fitness, 'history': self.history } ``` ### 贪心算法与构造式启发式 贪心算法简单快速，适合作为基准解： ```python import numpy as np from typing import List, Callable, Tuple, Dict class GreedyOptimizer: """ 贪心优化器 核心思想：每步选择局部最优 """ def __init__(self, candidates: List, evaluate_func: Callable): """ Args: candidates: 候选方案列表 evaluate_func: 评估函数，返回目标值 """ self.candidates = candidates self.evaluate_func = evaluate_func self.selected = [] self.history = [] def greedy_add(self, n_select: int, constraint_func: Callable = None) -> Tuple[List, float]: """ 贪心添加策略 每次选择能带来最大边际收益的候选 Args: n_select: 选择数量 constraint_func: 约束函数，返回True表示可行 Returns: (选择的候选列表, 最终目标值) """ available = set(range(len(self.candidates))) self.selected = [] for _ in range(n_select): best_candidate = None best_value = -float('inf') # 尝试每个可用候选 for idx in list(available): # 检查约束 test_selection = self.selected + [idx] if constraint_func and not constraint_func(test_selection): continue # 评估 value = self.evaluate_func(test_selection) if value > best_value: best_value = value best_candidate = idx if best_candidate is None: break # 没有可行的候选 # 选择最好的 self.selected.append(best_candidate) available.remove(best_candidate) self.history.append(best_value) return self.selected, best_value def greedy_remove(self, initial_solution: List, n_remove: int) -> Tuple[List, float]: """ 贪心移除策略 从初始解开始，每次移除损失最小的 Args: initial_solution: 初始解（候选索引列表） n_remove: 移除数量 Returns: (剩余候选列表, 最终目标值) """ current = set(initial_solution) self.selected = list(current) for _ in range(n_remove): if len(current) <= 1: break worst_candidate = None min_loss = float('inf') initial_value = self.evaluate_func(list(current)) # 尝试移除每个候选 for idx in list(current): test_selection = current - {idx} value = self.evaluate_func(list(test_selection)) loss = initial_value - value if loss < min_loss: min_loss = loss worst_candidate = idx if worst_candidate is not None: current.remove(worst_candidate) self.selected = list(current) self.history.append(self.evaluate_func(self.selected)) return self.selected, self.evaluate_func(self.selected) def adaptive_greedy(self, n_select: int, constraint_func: Callable = None) -> Tuple[List, float]: """ 自适应贪心 结合添加和移除策略，改进解质量 Args: n_select: 目标选择数量 constraint_func: 约束函数 Returns: (选择的候选列表, 最终目标值) """ # 先用贪心添加 selected, _ = self.greedy_add(n_select, constraint_func) # 尝试局部搜索改进 improved = True while improved: improved = False best_swap = None best_value = self.evaluate_func(selected) # 尝试交换 selected_set = set(selected) available_set = set(range(len(self.candidates))) - selected_set for out_idx in selected: for in_idx in available_set: new_selection = [in_idx if x == out_idx else x for x in selected] if constraint_func and not constraint_func(new_selection): continue value = self.evaluate_func(new_selection) if value > best_value: best_value = value best_swap = (out_idx, in_idx) if best_swap: out_idx, in_idx = best_swap selected = [in_idx if x == out_idx else x for x in selected] improved = True return selected, best_value ``` --- ## 代码示例 ### 生态网络优化问题 ```python """ 生态网络优化：在预算约束下最大化连通性 """ import numpy as np from typing import List, Tuple, Dict, Set import networkx as nx class EcologicalNetworkOptimizer: """ 生态网络优化器 目标：选择廊道建设方案，在预算约束下最大化生态连通性 """ def __init__(self, sources: List[Dict], corridor_candidates: List[Dict], budget: float): """ Args: sources: 生态源地列表 [{'id': i, 'pos': (x, y), 'quality': q}, ...] corridor_candidates: 候选廊道列表 [{'from': i, 'to': j, 'cost': c, 'quality': q}, ...] budget: 总预算 """ self.sources = sources self.candidates = corridor_candidates self.budget = budget self.n_sources = len(sources) self.n_corridors = len(corridor_candidates) # 构建源地图 self.source_map = {s['id']: i for i, s in enumerate(sources)} def build_graph(self, selected_corridors: List[int]) -> nx.Graph: """ 根据选中的廊道构建图 Args: selected_corridors: 选中的廊道索引列表 Returns: NetworkX图 """ G = nx.Graph() # 添加节点（源地） for source in self.sources: G.add_node( source['id'], pos=source['pos'], quality=source.get('quality', 1.0) ) # 添加边（廊道） for idx in selected_corridors: corridor = self.candidates[idx] G.add_edge( corridor['from'], corridor['to'], weight=corridor.get('quality', 1.0), cost=corridor['cost'], length=corridor.get('length', 1) ) return G def evaluate_connectivity(self, selected: List[int]) -> float: """ 评估连通性（目标函数） 综合考虑： 1. 连通源地数量 2. 最大连通分量大小 3. 网络平均最短路径 Args: selected: 选中的廊道索引列表 Returns: 连通性得分（越高越好） """ if not selected: return 0.0 G = self.build_graph(selected) if G.number_of_nodes() == 0: return 0.0 score = 0.0 # 1. 连通源地数 connected_sources = len([n for n in G.nodes() if G.degree(n) > 0]) score += connected_sources / self.n_sources * 0.4 # 2. 最大连通分量 if G.number_of_edges() > 0: largest_cc = max(len(cc) for cc in nx.connected_components(G)) score += largest_cc / self.n_sources * 0.4 else: score += 0.0 # 3. 网络效率（仅当图连通时） if nx.is_connected(G): # 使用平均最短路径的倒数（越短越好） avg_path = nx.average_shortest_path_length(G, weight='weight') efficiency = 1.0 / (1.0 + avg_path) score += efficiency * 0.2 return score def check_budget(self, selected: List[int]) -> bool: """检查是否满足预算约束""" total_cost = sum(self.candidates[i]['cost'] for i in selected) return total_cost <= self.budget def greedy_solve(self) -> Tuple[List[int], float]: """ 贪心算法求解 Returns: (选中的廊道索引列表, 得分) """ # 按性价比排序 candidates_with_idx = [ (i, c['cost'], c.get('quality', 1.0) / max(c['cost'], 1)) for i, c in enumerate(self.candidates) ] candidates_with_idx.sort(key=lambda x: -x[2]) # 按性价比降序 selected = [] remaining_budget = self.budget for idx, cost, _ in candidates_with_idx: if cost <= remaining_budget: selected.append(idx) remaining_budget -= cost # 评估 score = self.evaluate_connectivity(selected) return selected, score def genetic_solve(self, population_size: int = 100, n_generations: int = 200) -> Tuple[List[int], float]: """ 遗传算法求解 Returns: (选中的廊道索引列表, 得分) """ def objective_func(individual): """目标函数（最小化，所以取负）""" if not self.check_budget(individual): return 1e6 # 惩罚不可行解 return -self.evaluate_connectivity(individual) ga = GeneticAlgorithm( objective_func=objective_func, n_variables=self.n_corridors, variable_type='binary', population_size=population_size, mutation_rate=0.02, crossover_rate=0.8, elite_size=5 ) result = ga.evolve(n_generations) selected = [i for i, val in enumerate(result['best_solution']) if val == 1] score = self.evaluate_connectivity(selected) return selected, score def simulated_annealing_solve(self, max_iterations: int = 10000) -> Tuple[List[int], float]: """ 模拟退火求解 Returns: (选中的廊道索引列表, 得分) """ def objective_func(individual): """目标函数（最小化）""" if not self.check_budget([i for i, val in enumerate(individual) if val == 1]): return 1e6 return -self.evaluate_connectivity([i for i, val in enumerate(individual) if val == 1]) sa = SimulatedAnnealing( objective_func=objective_func, n_variables=self.n_corridors, variable_type='binary', initial_temp=100, cooling_rate=0.995, min_temp=0.01 ) result = sa.optimize(max_iterations) selected = [i for i, val in enumerate(result['best_solution']) if val == 1] score = self.evaluate_connectivity(selected) return selected, score def compare_methods(self) -> Dict: """ 比较不同求解方法 Returns: 各方法的结果 """ results = {} print("Running Greedy...") greedy_selected, greedy_score = self.greedy_solve() results['greedy'] = { 'selected': greedy_selected, 'score': greedy_score, 'cost': sum(self.candidates[i]['cost'] for i in greedy_selected) } print(f" Greedy: {len(greedy_selected)} corridors, score={greedy_score:.3f}") print("Running Genetic Algorithm...") ga_selected, ga_score = self.genetic_solve(population_size=50, n_generations=100) results['genetic'] = { 'selected': ga_selected, 'score': ga_score, 'cost': sum(self.candidates[i]['cost'] for i in ga_selected) } print(f" GA: {len(ga_selected)} corridors, score={ga_score:.3f}") print("Running Simulated Annealing...") sa_selected, sa_score = self.simulated_annealing_solve(max_iterations=5000) results['simulated_annealing'] = { 'selected': sa_selected, 'score': sa_score, 'cost': sum(self.candidates[i]['cost'] for i in sa_selected) } print(f" SA: {len(sa_selected)} corridors, score={sa_score:.3f}") return results # 使用示例 def example_network_optimization(): """生态网络优化示例""" # 创建生态源地 np.random.seed(42) n_sources = 10 sources = [ { 'id': i, 'pos': (np.random.uniform(0, 100), np.random.uniform(0, 100)), 'quality': np.random.uniform(0.5, 1.0) } for i in range(n_sources) ] # 创建候选廊道（所有源地对之间的连线） corridor_candidates = [] for i in range(n_sources): for j in range(i + 1, n_sources): pos_i = sources[i]['pos'] pos_j = sources[j]['pos'] length = np.sqrt((pos_i[0] - pos_j[0])**2 + (pos_i[1] - pos_j[1])**2) corridor_candidates.append({ 'from': sources[i]['id'], 'to': sources[j]['id'], 'cost': length * 10, # 成本与距离成正比 'quality': (sources[i]['quality'] + sources[j]['quality']) / 2, 'length': length }) # 设置预算 total_cost_all = sum(c['cost'] for c in corridor_candidates) budget = total_cost_all * 0.3 # 预算为总成本的30% # 创建优化器 optimizer = EcologicalNetworkOptimizer(sources, corridor_candidates, budget) # 比较方法 results = optimizer.compare_methods() # 找出最佳方法 best_method = max(results.keys(), key=lambda k: results[k]['score']) print(f"\nBest method: {best_method}") print(f"Best score: {results[best_method]['score']:.3f}") print(f"Corridors selected: {len(results[best_method]['selected'])}") print(f"Budget used: {results[best_method]['cost']:.1f} / {budget:.1f}") return optimizer, results if __name__ == "__main__": example_network_optimization() ``` --- ## 案例分析 ### ENAgent中的保护区优化 ENAgent使用空间优化算法设计最优的保护区网络： ```python class ReserveDesignOptimizer: """ 保护区设计优化器 基于 Marxan 思想：用最小成本实现保护目标 """ def __init__(self, planning_units: np.ndarray, features: Dict[str, np.ndarray], cost_surface: np.ndarray, targets: Dict[str, float]): """ Args: planning_units: 规划单元（可以是栅格） features: 各生态特征的分布 {'species': raster} cost_surface: 每个单元的保护成本 targets: 各特征的保护目标 {'species': proportion} """ self.planning_units = planning_units self.features = features self.cost_surface = cost_surface self.targets = targets self.n_units = planning_units.size # 计算每个特征的现有数量 self.feature_amounts = { name: (raster > 0).sum() for name, raster in features.items() } def objective_function(self, solution: np.ndarray) -> float: """ 目标函数：最小化成本 + 惩罚未达目标 Args: solution: 二进制解向量 Returns: 目标值（越小越好） """ # 成本 cost = (solution * self.cost_surface).sum() # 惩罚项 penalty = 0 penalty_factor = cost.sum() * 2 # 惩罚系数 for feature_name, raster in self.features.items(): # 计算被保护的特征量 protected_amount = (solution * (raster > 0)).sum() target_amount = self.feature_amounts[feature_name] * self.targets[feature_name] if protected_amount < target_amount: # 未达目标的惩罚 shortfall = target_amount - protected_amount penalty += shortfall * penalty_factor / self.feature_amounts[feature_name] # 边界长度惩罚（促进紧凑性） boundary_penalty = self._compute_boundary_length(solution) * 0.1 return cost + penalty + boundary_penalty def _compute_boundary_length(self, solution: np.ndarray) -> float: """计算边界长度（促进紧凑性）""" solution_2d = solution.reshape(self.planning_units.shape) # 计算边界 boundary = 0 for i in range(solution_2d.shape[0]): for j in range(solution_2d.shape[1]): if solution_2d[i, j] == 1: # 检查4邻域 for di, dj in [(-1, 0), (1, 0), (0, -1), (0, 1)]: ni, nj = i + di, j + dj if 0 <= ni < solution_2d.shape[0] and 0 <= nj < solution_2d.shape[1]: if solution_2d[ni, nj] == 0: boundary += 1 return boundary def iterative_improvement(self, initial_solution: np.ndarray, max_iterations: int = 1000) -> np.ndarray: """ 迭代改进算法 Args: initial_solution: 初始解 max_iterations: 最大迭代次数 Returns: 优化后的解 """ current = initial_solution.copy() current_value = self.objective_function(current) for iteration in range(max_iterations): improved = False # 尝试添加 for i in np.random.permutation(self.n_units): if current[i] == 0: test = current.copy() test[i] = 1 test_value = self.objective_function(test) if test_value < current_value: current = test current_value = test_value improved = True break if not improved: # 尝试移除 for i in np.random.permutation(self.n_units): if current[i] == 1: test = current.copy() test[i] = 0 test_value = self.objective_function(test) if test_value < current_value: current = test current_value = test_value improved = True break if not improved: break return current def solve(self, method: str = 'greedy') -> Tuple[np.ndarray, Dict]: """ 求解保护区设计问题 Args: method: 'greedy', 'iterative', 或 'simulated_annealing' Returns: (解, 结果信息) """ if method == 'greedy': return self._greedy_solve() elif method == 'iterative': return self._iterative_solve() elif method == 'simulated_annealing': return self._sa_solve() else: raise ValueError(f"Unknown method: {method}") def _greedy_solve(self) -> Tuple[np.ndarray, Dict]: """贪心求解""" # 按性价比排序 benefit_cost_ratio = np.zeros(self.n_units) for feature_name, raster in self.features.items(): feature_present = (raster > 0).flatten() benefit_cost_ratio += feature_present * self.cost_surface.flatten() benefit_cost_ratio = np.where( benefit_cost_ratio > 0, 1.0 / benefit_cost_ratio, 0 ) # 按性价比贪心选择 order = np.argsort(-benefit_cost_ratio) solution = np.zeros(self.n_units, dtype=int) current_cost = 0 max_cost = self.cost_surface.sum() * 0.3 # 预算约束 for idx in order: if benefit_cost_ratio[idx] > 0: test_cost = current_cost + self.cost_surface.flatten()[idx] if test_cost <= max_cost: solution[idx] = 1 current_cost = test_cost # 迭代改进 solution = self.iterative_improvement(solution) return solution, { 'cost': current_cost, 'objective': self.objective_function(solution), 'area_selected': solution.sum() } def _iterative_solve(self) -> Tuple[np.ndarray, Dict]: """迭代改进求解""" # 从贪心解开始 initial, _ = self._greedy_solve() solution = self.iterative_improvement(initial, max_iterations=1000) return solution, { 'cost': (solution * self.cost_surface.flatten()).sum(), 'objective': self.objective_function(solution), 'area_selected': solution.sum() } def _sa_solve(self) -> Tuple[np.ndarray, Dict]: """模拟退火求解""" def obj_func(x): return self.objective_function(x) sa = SimulatedAnnealing( objective_func=obj_func, n_variables=self.n_units, variable_type='binary', initial_temp=1000, cooling_rate=0.99, min_temp=0.1 ) result = sa.optimize(max_iterations=10000) return result['best_solution'], { 'cost': (result['best_solution'] * self.cost_surface.flatten()).sum(), 'objective': result['best_fitness'], 'area_selected': result['best_solution'].sum(), 'history': result['history'] } ``` --- ## 反思与延伸 ### 思考问题 1. **局部最优 vs 全局最优**：在空间优化中，局部最优解是否一定不可接受？ 2. **计算效率**：当问题规模达到百万级别时，如何平衡解质量和计算时间？ 3. **多目标权衡**：如何处理生态保护与经济发展的冲突？ 4. **不确定性**：数据不确定性如何在优化中考虑？ 5. **动态优化**：当环境条件变化时，如何更新优化解？ ### 延伸阅读 - **"Metaheuristics in Spatial Optimization"** - 空间优化综述 - **"Optimization Methods in GIS"** - GIS中的优化方法 - Marxan documentation - 保护区设计经典工具 - **"Integer Programming"** (Wolsey) - 整数规划理论 --- ## 关键要点 1. **空间优化 = 目标 + 约束 + 变量**：清晰的数学建模是成功的关键 2. **没有万能算法**：不同问题需要不同的求解策略 3. **精确算法适用于小规模**：大规模问题必须用启发式 4. **元启发式需要参数调优**：遗传算法、模拟退火等需要仔细设置参数 5. **解的稳健性很重要**：敏感性分析验证优化结果的可靠性