2026_DesignAI/officefile/md/supplements/02-spatial-intelligence/02.4-spatial-optimization.md

# 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. **解的稳健性很重要**：敏感性分析验证优化结果的可靠性