pengxiao
a90f7adfa1
将 Markdown 源文件移入 md/，LaTeX 工作目录保留在 latex/， Word 导出移入 word/；删除临时脚本、调试截图和空 stub。 Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
45 KiB
45 KiB
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)
遗传算法模拟自然进化过程，是空间优化中最常用的元启发式方法：
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)
模拟退火模拟金属冷却过程，能跳出局部最优：
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
}
贪心算法与构造式启发式
贪心算法简单快速，适合作为基准解：
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
代码示例
生态网络优化问题
"""
生态网络优化：在预算约束下最大化连通性
"""
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使用空间优化算法设计最优的保护区网络：
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']
}
反思与延伸
思考问题
-
局部最优 vs 全局最优：在空间优化中，局部最优解是否一定不可接受？
-
计算效率：当问题规模达到百万级别时，如何平衡解质量和计算时间？
-
多目标权衡：如何处理生态保护与经济发展的冲突？
-
不确定性：数据不确定性如何在优化中考虑？
-
动态优化：当环境条件变化时，如何更新优化解？
延伸阅读
- "Metaheuristics in Spatial Optimization" - 空间优化综述
- "Optimization Methods in GIS" - GIS中的优化方法
- Marxan documentation - 保护区设计经典工具
- "Integer Programming" (Wolsey) - 整数规划理论
关键要点
-
空间优化 = 目标 + 约束 + 变量：清晰的数学建模是成功的关键
-
没有万能算法：不同问题需要不同的求解策略
-
精确算法适用于小规模：大规模问题必须用启发式
-
元启发式需要参数调优：遗传算法、模拟退火等需要仔细设置参数
-
解的稳健性很重要：敏感性分析验证优化结果的可靠性