2407054130黄联鑫B题 Jupyter Notebook PyCharm Community Edition

[乂忆]

土壤普查路径和任务规划

摘   要

本文针对土壤普查路径和任务规划问题,建立了一系列数学模型,并使用多种优化算法和数据处理技术进行求解,以期在有限的时间和资源条件下,高效完成大量采样点的采样任务。

在问题1的求解过程中，根据“工作组在所有点位之间可以按照直线通行”的假设采用半正矢公式计算两地间距离，然后采用采用动态规划的Held-Karp算法求解八个点间最短哈密顿路径。通过图表来直观呈现结果。

在问题2的求解过程中，需要依据就近原则将222个点以8为基数进行划分为27个簇，采用k-means聚类算法进行划分。每个簇代表一天的工作任务，依据问题1的算法计算最优路径与时间。通过图表来直观呈现结果。

在问题3的求解过程中，通过构建条件判断来求解第一问是否出现超过时限，根据方差公式计算是否存在工作时间不均衡现象，问题二通过整数规划的模型来进行建模， 以每天工作的时间作为约束条件，然后使用分支定界法对其进行求解得到均衡化方案。通过图表来直观呈现结果。

在问题4的求解过程中，利用第三方地图服务提供的路径规划API,获取采样点之间的实际行驶距离和时间数据,替换问题一中所使用的距离数据集，得到考虑实际路况的最优采样路径。通过图表来直观呈现结果。

关键词：

土壤普查路径和任务规划 k-means聚类算法 最短哈密顿路径 整数规划 地图 API Held-Karp算法

 

 

土壤普查路径和任务规划

 

• ．问题重述

土壤普查是确保国家粮食安全的基础工作之一。研究土壤普查路径和任务规划具有一定现实意义。根据题目要求，需要解决以下几个问题：

（1）：用经纬度距离替代两点距离的方式，以速度20千米/小时计算赶路时间，求八个工作点间的最短哈密顿路径和最短工作时间。

（2）：采用就近原则，以每天工作8个点位为基准，原222个点仅用前216点划分为27个簇，用经纬度距离替代两点距离的方式，以速度20千米/小时计算赶路时间。求每一天的点位最优路径，及相应的工作时间，并对所有的最优路径工作时间求取最大值和最小值。

（3）：计算最优路径工作时间是否超过8.5小时，并计算每天方差的大小变化幅度，判断是否工作超时，方案是否均衡。取消每天8个点位数的限制，在工作时限内求均衡化最优路径方案。

（4）：用地图API计算专车通过8个点位的时间，求最优路径和最短工作时间。

（二）．问题的分析

问题一要求用经纬度距离替代两点距离的方式，以速度20千米/小时计算赶路时间，求八个工作点间的最短哈密顿路径和最短工作时间。对于此问题，我们认为关键是计算两点间距离、求最短哈密顿路径，为了解决本问题，我们拟采用半正矢公式求8个点间的距离，用动态规划的Held-Karp算法求最短哈密顿路径我们认为该问题本质是求最短哈密顿路径。

问题二要求采用就近原则，以每天工作8个点位为基准，将222个点划分为27个簇，求每一天的点位最优路径，及相应的工作时间，并对所有的最优路径工作时间求取最大值和最小值。对于此问题，我们认为关键是点位的划分为了解决本问题，我们拟采用k-means聚类算法进行划分簇。往后的问题是对上一题的简单重复，用先前算法可以解。

问题三要求算最优路径工作时间是否超过8.5小时，判断是否工作超时，方案是否均衡。取消每天8个点位数的限制，在工作时限内求均衡化最优路径方案。对于此问题，我们认为关键是取消点位数限制后，求工作时限内均衡化最优路径方案为了解决本问题，我们拟采用分支定界法来获得最优路径。

问题四要求用地图API计算专车通过8个点位的时间，求最优路径和最短工作时间对于此问题，我们认为关键是用地图API计算通过点位的时间为了解决本问题，我们拟采用第三方地图服务提供的路径规划API求距离，用问题一中算法求解后续问题。

（三）．模型假设

假设一：问题一、二、三，两点间路径为假设为直线距离，速度恒定，工作时间恒定，不考虑路况、天气，工作效率等细微因素。

假设二：问题四中，由于采用专车通行的方案，直线假设不能使用，需要用第三方地图软件提供API进行路径时间测算，所求时间均为最短时间，不考虑交通成本。

假设三：在问题二的模型建立过程中，论文假设k-means 聚类算法可以成功划分所有采样点,每个簇代表一天的工作任务，不考虑偏离点的问题，将各点按就近原则分成 27 个簇

假设四：问题3,论文假设可以建立整数规划模型来描述问题,其中决策变量为二进制变量,表示每个采样点是否被分配到某一天;目标函数为去除8点限制后，最小化完成所有采样任务的最长工作时间,约束条件包括每个采样点只能被分配到一天、 每天的工作时间不超过 8.5 小时等;在求解该整数规划模型时,论文假设可以使用分支定界法等经典算法来获得全局最优解。

（四）．符号及变量说明

Lng1 Lng2表示经度值

Lat1 Lat2表示维度值

 V={0,1,2,…,n−1}，表示 n 个需要访问的点位。

S：表示已访问的城市集合，用位掩码表示

dp(S,i)：表示访问子集 S 中的城市且最后一个到达城市 i 的最短路径长度。

d(i,j)：表示城市 i 到城市 j 的移动时间（或距离）。

是第K个簇的点集。

是属于的数据点。

是第K个簇的质心。

 表示数据点与簇中心之间的欧氏距离平方。

∑求和符号
（五）．模型的建立与求解

问题一：最短哈密顿路径

1、建模思路分析

问题需要求一条最短哈密顿路径，预计采用动态规划的Held-Karp算法算法，对于路径采用haversine公式，由于该问题与后续问题有关联，方法可重复利用，预计编写python代码编写。

1. 两点间距离测算

两点间距离测算需考虑地球形状，采用半正矢公式进行基于经纬度的距离测算，半正矢公式如下。

其中：

1. Lng1 Lat1 表示A点经纬度，Lng2 Lat2 表示B点经纬度

2. a=Lat1 - Lat2 为两点纬度之差 b=Lng1 - Lng2 为两点经度之差

3. 6378.137为地球赤道半径，单位为千米

下图为两点间距离，单位为千米，A-G分别代表1-8，点位依次为31、44、61、83、100、115、147、158。

由此可得通行时间如下，单位为分钟。

1. 求最短哈密顿路径

预采用Held-Karp算法递归地解决子问题，并存储中间结果以避免重复计算。

公式如下

对于每个单点集合 S={i}，路径长度为0（无移动）：

dp({i},i)=0 ∀i∈V

    对于每个子集 S（∣S∣≥2）和每个城市 i∈S，通过遍历所有可能的上一城市 j∈S∖{i}，计算最短路径：

dp(S,i)=min（j∈Sj≠i）[dp(S∖{i},j)+d(j,i)]

访问所有城市后的最短开放路径长度（无需返回起点）：

最短移动时间= min（i∈V）dp(V,i)

总工作时间包括所有点位的工作时间和移动时间之和：

总时间= ∑（i∈V）t work(i)+最短移动时间

其中 t work(i) 为点位 i 的工作时间。

由此可得哈密顿路径如下

哈密顿路径为115、31、100、147、83、158、61、44

所求最短时间为400.10分钟，约6.67小时，处于理想环境，若考虑其他实际因素，最短时间可能延长。

通过Held-Karp算法，我们成功找到一条哈密顿路径，缩短了行路时间，降低了成本，该模型基于经纬度距离建立，若要提高精确度，需要采用地图API算法与更加精确的经纬度距离算法。

问题二：聚类分组

1. 建模思路分析

该问题要求将附件中222个点以8个点为基准分为27个簇，给出每一天的点位最优路径，及相应的工作时间，并对所有的最优路径工作时间求取最大值和最小值。

该问题可以分为两部分，一：将222个点位划分，二：求每一天8个点位的最优路径与最小时间。

对于第一部分，可以采用k-means聚类算法，将相近点位划分为一个簇，以缩短工作时间。

对于第二部分，其可以通过修改点位化归为问题一，计算每个簇中的路径与时间，并求最大值和最小值，依旧是求最小哈密顿路径。

2、k-means聚类模型建立

将222个点位进行分组，可以采用k-means聚类算法，其步骤如下：

初始化：随机选择K个数据点作为初始簇质心。

分配：计算每个点位与数据点位的距离，将每个数据点分配给最近的簇质心。

更新：重新计算每个簇的质心，通常是簇内所有数据点的均值。

迭代：重复分配和更新步骤，直到满足终止条件，如簇质心的变化小于某个阈值或达到预设的迭代次数。

K-means 的数学公式

K-means 的目标是最小化簇内平方误差和，即每个点到其所属簇中心的距离的平方和，公式如下：

其中：

 K是簇的数量。

是第K个簇的点集。

是属于的数据点。

是第K个簇的质心。

 表示数据点与簇中心之间的欧氏距离平方。

可以直接计算

Python中使用约束聚类库k-means-constrained进行聚类

1. 结果呈现

第一组结果如下

第二组结果如下

第三组结果如下

第四组结果如下

第五组结果如下

第六组结果如下

第七组结果如下

第八组结果如下

第九组结果如下

第十组结果如下

第十一组结果如下

第十二组结果如下

第十三组结果如下

第十四组结果如下

第十五组结果如下

第十六组结果如下

第十七组结果如下

第十八组结果如下

第十九组结果如下

第二十组结果如下

第二一组结果如下

第二二组结果如下

第二三组结果如下

第二四组结果如下

第二五组结果如下

第二六组结果如下

第二七组结果如下

总结路径，时间，如下

最大时间为第20组，时间为450.46分钟，约为7.51小时

最小时间为第10组，时间为386.15分钟，约为6.44小时

    如上为有拘束的k-means聚类与Held-Karp算法所解出的任务划分和路径优化问题。 虽然每个工作日的具体情况有所不同,但通过合理的任务分配和路径规划,我们最终得到了一个较为均衡的工作方案,满足了实际工作的需求。

问题三：平衡工作时间

1. 问题分析

问题 3 要求在每天工作时间不超过 8.5 小时的条件下,给出一个更加均衡的工作

安排方案,使得尽快完成所有采样任务。

这个问题的关键在于如何合理地分配222个采样点的工作任务,在满足每天工作时间限制的情况下,尽量平衡各天的工作量。

我们可以考虑使用一种优化算法,如整数规划或混合整数规划,来建立数学模型并求解这个问题。整数规划模型可以很好地描述问题中的离散决策变量和约束条件。

1. 数据分析

根据第二题中6.44-7.51小时可得

数据可得问题二所给的方案不存在工作组每天工作时间超过限时现象

工作时间最大差值为1.07小时

工作时间较为均衡

1. 优化分配

根据各个采样点的完成工作所需时间,构建整数规划模型。

为求解这个整数规划模型,可以使用一些常见的优化算法,如分支定界法、切平面法或近似算法等。

以下是使用分支定界法求解问题 3 的思路:

(1) 初始化:设置问题的上下界,即完成所有采样任务的最短和最长工作时间。

(2) 分支:选择一个未固定的决策变量�𝑖,将问题分为两个子问题:�𝑖 = 0 和�𝑖 = 1。

(3) 定界:求解每个子问题的下界,并与当前最优解进行比较,剪掉不可能达到最优解的子问题。

(4) 选择:选择一个未固定的决策变量,重复步骤(2)和(3),直到所有变量都被固定。

(5) 更新:更新当前最优解。

(6) 回溯:如果当前分支无法得到更优解,则回溯到上一层分支,选择另一个分支继续探索。

(7) 终止:当所有分支都被探索完毕,或达到预设的终止条件(如计算时间或迭代次数),算法终止。输出最优解:输出完成所有采样任务的最长工作时间,以及每个采样点被分配到的工作日。

在分支定界法中,我们需要解决以下几个关键问题:

如何确定问题的上下界:

可以通过简单的启发式方法,如将所有点位平均分配到 k 天,或者将所有点位按照最短路径依次分配,来确定初始的上下界。

如何选择分支变量:

可以采用启发式策略,如选择当前最长工作时间对应的采样点,或者选择当前未分配

的采样点中工作时间最长的那个。

如何计算子问题的下界:

可以使用松弛问题或者启发式方法,如将当前已分配的点位固定,然后求解剩余点位

的最优分配。

4、算法分析

分支定界法是求解整数规划问题的经典算法之一。 它通过不断地对问题进行分

支和定界,来缩小可行解空间,最终找到全局最优解。

问题四：专车通行方案

1. 问题分析

该问题需要将第一问中的直线假设去掉，改为用地图API计算通行时间，后求最短哈密顿路径的方法不变。

1. 求通行时间

通过调用百度API获得数据如下

单位为分钟

将数据带入问题一中代码可以获得哈密顿路径与时间

获得结果如下

最优路径顺序：[31, 115, 100, 147, 83, 158, 61, 44]

最短总时间：421.00 分钟

最短总时间：7.02 小时

总结：

    使用地图API为理想模型，不考虑其他干扰性因素，实际情况需实际测量，存在范围差值。

（六）模型的检验、灵敏度分析

对于该问题，目前算法以可以大致解决，若需更精确的模型，需要更改精确的数据测量工具，可以使用地图API进行距离测算。

在目前，由于haversine公式存在一定误差，故时间不一定准确，在灵敏度上，随地球半径的变化较大，需精确。

（七）模型的优缺点分析，模型的推广

优点：可以很好的解决题设问题，便于推广到更多的哈密顿路径问题

缺点：存在理想假设，如需投入使用，需要考虑干扰性因素，建议实际测量通行时间为准。

（八）参考文献

（九）附录

问题一代码：

import math

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

from itertools import combinations

 

pro1=pd.read_excel('pro1.xlsx')

pro1

 

 

 

S=np.zeros((8,8))

for x in range(0,8):

    for y in range(0,8):

        a=math.radians(pro1['WD'][x]-pro1['WD'][y])

        b=math.radians(pro1['JD'][x]-pro1['JD'][y])

        S[x][y]=2*6378.137*math.asin(math.sqrt((math.sin(a/2))**2+((math.sin(b/2))**2)*math.cos(math.radians(pro1['WD'][x]))*math.cos(math.radians(pro1['WD'][y]))))

 

n = 8

work_total = sum(pro1['TIME'])

 

# 构建时间矩阵

time_matrix = [[0]*n for _ in range(n)]

for i in range(n):

    for j in range(n):

        if i != j:

            time_matrix[i][j] = S[i][j]  # 小时

 

# 动态规划求解TSP

INF = float('inf')

dp = [[INF]*n for _ in range(1 << n)]

predecessor = [[None]*n for _ in range(1 << n)]

 

for u in range(n):

    dp[1 << u][u] = 0

 

masks = sorted([(bin(mask).count('1'), mask) for mask in range(1 << n)], key=lambda x: x[0])

 

for count, mask in masks:

    if count == 0:

        continue

    for u in range(n):

        if not (mask & (1 << u)):

            continue

        if dp[mask][u] == INF:

            continue

        for v in range(n):

            if mask & (1 << v):

                continue

            new_mask = mask | (1 << v)

            new_cost = dp[mask][u] + time_matrix[u][v]

            if new_cost < dp[new_mask][v]:

                dp[new_mask][v] = new_cost

                predecessor[new_mask][v] = u

 

full_mask = (1 << n) - 1

min_time = min(dp[full_mask][v] for v in range(n))

best_v = dp[full_mask].index(min_time)

 

 

# 回溯路径

path = []

current_mask, current_v = full_mask, best_v

while current_mask != (1 << current_v):

    path.append(current_v)

    prev_v = predecessor[current_mask][current_v]

    current_mask ^= (1 << current_v)

    current_v = prev_v

path.append(current_v)

path = path[::-1]

 

o=-1

mint=0

for i in path:

    if(o==-1):

        o=i;

    else:

        mint=S[o][i]+mint

 

# 结果输出

optimal_order = [pro1['ID'][i] for i in path]

total_time = work_total + mint

time=total_time

print(f"最优路径顺序：{optimal_order}")

print(f"最短总时间：{time:.2f} 分钟")

time=time/60

print(f"最短总时间：{time:.2f} 小时")

plt.figure(figsize=(5, 4))

 

# 绘制城市位置

plt.scatter(pro1['JD'], pro1['WD'], c='red', s=100,

           edgecolors='black', zorder=3)

   

# 标注城市名称（假设有'City'列）

for i in range(len(pro1)):

    plt.text(pro1['JD'][i]+0.0005, pro1['WD'][i]+0.0005,

                pro1['ID'][i], fontsize=10, ha='left')

   

# 绘制最优路径

ordered_jd = [pro1['JD'][i] for i in path]

ordered_wd = [pro1['WD'][i] for i in path]

   

plt.plot(ordered_jd, ordered_wd, 'b--', alpha=0.6)

plt.quiver(ordered_jd[:-1], ordered_wd[:-1],

            np.diff(ordered_jd), np.diff(ordered_wd),

          scale_units='xy', angles='xy', scale=1,

          color='green', width=0.003, headwidth=8)

   

plt.grid(True, alpha=0.3)

plt.tight_layout()

plt.show()

问题二代码

import pandas as pd

from k_means_constrained import KMeansConstrained  # 注意库名不同

 

# 读取数据（假设坐标列名为'X'和'Y'）

df = pd.read_excel("dtd.xlsx")

 

# 约束聚类：27个簇，每个簇严格8个点

model = KMeansConstrained(

    n_clusters=27,

    size_min=8,

    size_max=8,

    random_state=42

)

 

# 执行聚类并保存结果

coordinates = df[['JD', 'WD']].values

model.fit(coordinates)

df['Cluster'] = model.labels_

df.to_excel("clustered_output.xlsx", index=False)

问题四地图API调用代码

import pandas as pd

import requests

import time

import numpy as np

from tqdm import tqdm  # 进度条显示

 

# 配置参数

INPUT_FILE = "pro1.xlsx"

OUTPUT_MATRIX_FILE = "time_matrix.xlsx"

AK = "     "  # 替换真实AK

COORD_TYPE = "wgs84"  # 坐标系类型

API_URL = "http://api.map.baidu.com/direction/v2/driving"

REQUEST_INTERVAL = 0.5  # 请求间隔(秒)

 

def get_driving_time(origin, destination, ak):

    """获取两点间行车时间"""

    params = {

        "origin": f"{origin[0]},{origin[1]}",

        "destination": f"{destination[0]},{destination[1]}",

        "ak": ak,

        "coord_type": COORD_TYPE,

        "tactics": 11  # 常规路线

    }

   

    try:

        response = requests.get(API_URL, params=params, timeout=10)

        result = response.json()

       

        if result.get('status') == 0:

            return result['result']['routes'][0]['duration']  # 返回秒数

        else:

            return f"Error {result.get('status')}"

    except Exception as e:

        return f"Exception: {str(e)}"

 

# 读取地点数据

df = pd.read_excel(INPUT_FILE)

locations = list(zip(df['WD'], df['JD']))  # 假设列名为纬度和经度

n = len(locations)

 

# 初始化时间矩阵

time_matrix = np.zeros((n, n), dtype=object)

location_names = df['ID'].tolist()  # 假设有地点名称列

 

# 生成进度条

pbar = tqdm(total=n*n, desc="生成时间矩阵")

 

# 遍历所有地点组合

for i in range(n):

    for j in range(n):

        if i == j:

            time_matrix[i][j] = 0  # 相同地点设为0

        else:

            time_matrix[i][j] = get_driving_time(locations[i], locations[j], AK)

            time.sleep(REQUEST_INTERVAL)

        pbar.update(1)

 

pbar.close()

 

# 转换为DataFrame保存

matrix_df = pd.DataFrame(time_matrix,

                        index=location_names,

                        columns=location_names)

 

matrix_df.to_excel(OUTPUT_MATRIX_FILE)

print(f"时间矩阵已生成至 {OUTPUT_MATRIX_FILE}")

[釉色清风]

不错