关于1999C煤矸石堆积的研读报告

姜智浩 Lv5

题目背景

此竞赛题目解决的是煤矸石堆积的问题 在面对矸石山不断堆积的情况下 需要制定合理的年度征地计划和给出在不同出矸率的情况下处理矸石的最低费用

模型假设

为了简化问题 此文章做了如下假设
photo

我认为这些假设是比较合理的 通过查询相关资料 其假设的数据符合实际生活中的数据

建模准备

在建模前他们首先给出了轨道长度与体积的关系、矸石山占地面积与体积的关系、占地面积与安息角的关系 在阅读的时候 我一开始不理解为什么要列出这些公式 在后续的阅读中我了解到了 轨道长度决定了运矸车运输的距离 从而影响能耗 成本 矸石堆积体积越大 轨道需要延伸得越长 建立轨道长度与体积的关系是为了后续计算运输成本 能量消耗等提供基础 占地面积决定率征地需求

征地费用是主要经济支出之一 因此需要知道给定体积下 占地多少土地 安息角决定了矸石自然堆积的最大坡度 坡度不同 相同体积下的占地面积就不同
我使用python对轨道长度与体积的关系、矸石山占地面积与体积的关系、占地面积与安息角的关系进行了复现 在假设体积为1000 安息角为35 轨道倾角为5时 得到轨道长度为61.3176 占地面积为35853292.5453 占地面积与安息角的关系如图所示:
photo

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
V = 1000  # 矸石山体积
beta_deg = 1 # 轨道倾角(度)
alpha_deg_range = np.linspace(1, 55, 500) # 安息角范围

beta = np.deg2rad(beta_deg)
alpha_rad = np.deg2rad(alpha_deg_range)

def compute_a(alpha, beta):
tan_alpha_sq = np.tan(alpha)**2
tan_beta_sq = np.tan(beta)**2
return tan_beta_sq / (tan_alpha_sq - tan_beta_sq)

def compute_a1(a):
sqrt_a = np.sqrt(a)
return 1 + (np.pi / 2 + np.arctan(sqrt_a)) * sqrt_a

def compute_S(V, alpha, beta):
a = compute_a(alpha, beta)
a1 = compute_a1(a)
numerator = 9 * a1 * V**2
denominator = np.sqrt(a) * np.tan(alpha)**2
S = (numerator / denominator)**(1/3)
return S

S_values = compute_S(V, alpha_rad, beta)

plt.figure(figsize=(10, 6))
plt.plot(alpha_deg_range, S_values, label=r'$S(\alpha)$', color='blue')
plt.xlabel(r'$\alpha$')
plt.ylabel('S')
plt.title(r'S VS $\alpha$')
plt.grid(True)
plt.show()

这与论文中的图标几乎一致 论文中的\alpha范围应该进行了归一化处理 我这里没有进行处理 有些许不同 但无关紧要

问题一

为了制定合理的年度征地计划 该论文建立了土地需求量公式 地价公式 电费公式 每年计划征地经费公式 征地规划模型
这些模型不是孤立存在的 而是环环相扣的 其最终目标优化土地使用 控制成本 合理安排资金流

土地需求量公式

在论文中 通过几何分析与经验拟合 得到矸石山占地面积:

1
S_i\ =\ 0.0015\left(\frac{9a_ii^2V_0^2}{\sqrt a{tan}^2\alpha}\right)^\frac{1}{3}

此公式中参数a和a1来源于对堆积体形状的几何分析
a受\alpha,β的影响 反映了轨道倾角\beta与安息角\alpha对堆积形态的影响 a1可以看成一个修正因子 用于描述此煤矸石底面积变化趋势 0.0015是一个系数 由于体积的单位是立方米 而占地面积要转换成亩 因此乘以0.0015用于系数转换
∆Si = Si - Si-1
这是一个显而易见的公式 表明了增加的变化量
因此土地需求量公式为
Si‘ = 1+b∆Si
通过python复现当\beta=25°时S1‘为10.663024317273257这与论文中的数值是一样的
不过我认为论文中的\beta=25°不妥当此数值设置得并不合理 可能要讨论不同仰角的变化 不过在后续的灵敏度分析中讨论了不同轨道仰角的效率和经费剩余

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
def compute_a(alpha_deg, beta_deg):
alpha = np.deg2rad(alpha_deg)
beta = np.deg2rad(beta_deg)
tan_alpha_sq = np.tan(alpha)**2
tan_beta_sq = np.tan(beta)**2
a = tan_beta_sq / (tan_alpha_sq - tan_beta_sq)
return a

def compute_a1(a):
sqrt_a = np.sqrt(a)
a1 = 1 + (np.pi / 2 + np.arctan(sqrt_a)) * sqrt_a
return a1

def compute_G0(d, e):
'''
求g0:每年产生的煤矸石质量
d:原煤年产量
e:出矸率
'''
g0 = (10000 * d * e) / (1 - e)
return g0

def compute_V0(g0, rho):
'''
求v0:每年产生的煤矸石体积
g0:每年产生的煤矸石质量
rho:煤矸石容量
'''
v0 = g0 / rho
return v0
def compute_Vi(i, v0):
'''
求vi:第i年产生的煤矸石体积
i:第i年
v0:每年产生的煤矸石体积
'''
vi = i * v0
return vi
def compute_Si(i, v0, alpha, beta):
'''
求Si:第i年煤矸石占地面积
i:第i年
v0:每年产生的煤矸石体积
alpha:安息角

a:计算参数a 根据compute_a可以计算得到
a1:计算参数a1 根据compute_a1可以计算得到
'''
a = compute_a(alpha, beta)
a1 = compute_a1(a)
tan_alpha = np.tan(np.deg2rad(alpha))
si = 0.0015 * ((9 * a1 * i**2 * v0**2) / (np.sqrt(a) * tan_alpha**2)) ** (1/3)
return si

def compute_delta_Si(i, v0, alpha, beta):
'''
求第i年末煤矸石占地面积的增量
'''
if i == 0:
return 0
else:
return compute_Si(i, v0, alpha, beta) - compute_Si(i - 1, v0, alpha, beta)

def compute_S_dot_i(i, v0, alpha, beta, b):
'''
求第i年土地需求量
'''
delta_Si = compute_delta_Si(i, v0, alpha, beta)
s_dot_i = (1 + b) * delta_Si
return s_dot_i

# 参数
beta = 25 # 轨道倾角
d = 300 # 原煤年产量
e = 0.07 # 出矸率
rho = 2 # 煤矸石容量
alpha = 55 # 安息角
i = 1 # 第 i 年
b = 0.1# 土地预留系数
g0 = compute_G0(d, e)
v0 = compute_V0(g0, rho)

result = compute_S_dot_i(i, v0, alpha, beta, b)
print(result)

def compute_C1i(i):
'''
求第i年的地价
'''
c0 = 8 #手动输入 根据假设6给定
lam = 0.1 #手动输入 根据假设6给定
c1i = c0 * (1 + lam) ** (i - 1)
return c1i

# 参数
i = 1 # 第 i 年
c1i = compute_C1i(i)
print(c1i)

征地规划模型

通过复现征地规划模型 我们得到了其结果 这个结果与论文中的一致 值得注意的是 在年数少于20时 其前几年得到的结果与论文中的不一致 经过排查得到在不同年数时资金时间价值累积效应 贷款安排策略 囤地策略会受到影响 本论文并未讨论不同年数对征地规划的影响 我认为此讨论是有必要的 为了增强论文的鲁棒性以及案例的迁移性 应该增加不同年份的征地策略 经过我的验证 当n=5时 只有前两年会选择贷款 并且每年的累计余额也不同 随着年数的增加 贷款会越来越激进 说明征地的短期与长期规划是截然不同的 在前期时利润是偏低的 但到了后期利润是巨大的
photo

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
import pulp

# 参数
n = 20
R1 = 0.05 # 自有资金利率
R2 = 0.05 # 贷款利率
c1 = [8, 8.8, 9.68, 10.648, 11.7128,
12.88, 14.17, 15,59, 17,15, 18.86,
20.75, 22.82, 25.11, 27.62, 30.38,
33.42, 36.76, 40.44, 44.48, 48.93] # 每年土地价格(万元/亩)

f1 = [97.11, 92.33, 90.31, 88.68, 87.26,
85.97, 84.78, 83.67, 82.61, 81.6,
80.62, 79.69, 78.78, 77.89, 77.03,
76.19, 75.37, 74.56, 73.77, 72.99] # 每年收入(万元)

S_given = [10.66, 6.26, 5.25, 4.69, 4.31,
4.03, 3.81, 3.63, 3.48, 3.36,
3.25, 3.15, 3.06, 2.99, 2.92,
2.85, 2.79, 2.74, 2.69, 2.64] # 每年实际围地面积

# 模型
model = pulp.LpProblem("1999C_Coal_Gangue_Model", pulp.LpMaximize)

# 决策变量
g = pulp.LpVariable.dicts("g", range(n + 1), lowBound=0) # 实际支出
u = pulp.LpVariable.dicts("u", range(n + 1), lowBound=0) # 贷款金额
w = pulp.LpVariable.dicts("w", range(n + 1), lowBound=0) # 征地费用
S_star = pulp.LpVariable.dicts("S_star", range(n + 1), lowBound=0) # 累计囤地
p = pulp.LpVariable.dicts("p", range(n + 1), lowBound=0) # 经费余额

# 初始约束
model += p[0] == 0, "初始经费"
model += S_star[0] == 0, "初始囤地"
model += u[0] == 0, "初始贷款"
model += u[n] == 0, "最后一年不贷款"
model += S_star[n] == 0, "最终囤地为0"

# 目标函数
model += p[n], "最终经费最大"

# 逐年约束
for i in range(1, n + 1):
# 征地费用计算
model += w[i] == g[i] - u[i - 1] * (1 + R2) + u[i], f"征地费用_{i}"

# 囤地数量递推
model += S_star[i] == S_star[i - 1] + w[i] * (1.0 / c1[i - 1]) - S_given[i - 1], f"囤地数量递推_{i}"

# 经费余额
model += p[i] == (p[i - 1] + f1[i - 1] - g[i]) * (1 + R1), f"经费累计余额_{i}"

# 支出限制
model += g[i] <= p[i - 1] + f1[i - 1], f"支出限制_{i}"

# 不能超额围地
for i in range(1, n):
model += S_star[i] <= sum(S_given[k] for k in range(i, n)), f"围地策略_{i}"

# 求解
model.solve(pulp.PULP_CBC_CMD(msg=0))

if pulp.LpStatus[model.status] != 'Optimal':
print("求解失败,状态:", pulp.LpStatus[model.status])
else:
print("求解成功,目标函数值(第20年末经费余额):", round(p[n].varValue, 2), "万元")
for i in range(1, n + 1):
print(f"\n第{i}年:")
print(f" 支出 g_{i} = {g[i].varValue:.2f}")
print(f" 贷款 u_{i} = {u[i].varValue:.2f}")
print(f" 囤地 S*_i = {S_star[i].varValue:.2f}")
print(f" 征地费用 w_{i} = {w[i].varValue:.2f}")
print(f" 经费余额 p_{i} = {p[i].varValue:.2f}")

问题二

我们同样复现了此模型 得到:
Optimal - objective value -35.832721
这个结果与论文中的一致 同样 我们发现 在不同的年数其最优解存在不同因此 我认为该论文添加上对不同年份的经费讨论

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
import pulp
import numpy as np

# 参数
d = 3000000 # 原煤年产量(吨)
rho = 2 # 煤矸石容量(吨/立方米)
alpha = 55 # 安息角(度)
beta = 25 # 轨道倾角(度)
n = 20 # 规划年限
R1 = 0.05 # 经费累计余额利率
R2 = 0.05 # 贷款利率

c1 = [8, 8.8, 9.68, 10.648, 11.7128,
12.88, 14.17, 15,59, 17,15, 18.86,
20.75, 22.82, 25.11, 27.62, 30.38,
33.42, 36.76, 40.44, 44.48, 48.93] # 每年土地价格(万元/亩)

f1 = [97.11, 92.33, 90.31, 88.68, 87.26,
85.97, 84.78, 83.67, 82.61, 81.6,
80.62, 79.69, 78.78, 77.89, 77.03,
76.19, 75.37, 74.56, 73.77, 72.99] # 每年收入(万元)

S_prime = [10.66, 6.26, 5.25, 4.69, 4.31,
4.03, 3.81, 3.63, 3.48, 3.36,
3.25, 3.15, 3.06, 2.99, 2.92,
2.85, 2.79, 2.74, 2.69, 2.64] # 每年实际围地面积

# 计算 a 和 a1
def compute_a(alpha_deg, beta_deg):
alpha = np.deg2rad(alpha_deg)
beta = np.deg2rad(beta_deg)
tan_alpha_sq = np.tan(alpha)**2
tan_beta_sq = np.tan(beta)**2
a = tan_beta_sq / (tan_alpha_sq - tan_beta_sq)
return a

def compute_a1(a):
sqrt_a = np.sqrt(a)
a1 = 1 + (np.pi / 2 + np.arctan(sqrt_a)) * sqrt_a
return a1

# 计算第 i 年的占地面积 Si
def compute_Si(i, v0, alpha):
a = compute_a(alpha, beta)
a1 = compute_a1(a)
si = 0.0015 * ((9 * a1 * i**2 * v0**2) / (np.sqrt(a) * np.tan(np.deg2rad(alpha))**2)) ** (1/3)
return si

def solve_model(v0, e):
# 计算第一年占地 S1
S1 = compute_Si(1, v0, alpha)

model = pulp.LpProblem("Minimum_Delta_q", pulp.LpMinimize)

# 决策变量
delta_q = pulp.LpVariable("delta_q", lowBound=None, upBound=None)
g_vars = pulp.LpVariable.dicts("g", range(1, n+1), lowBound=0)
u_vars = pulp.LpVariable.dicts("u", range(n+1), lowBound=0)
S_star = pulp.LpVariable.dicts("S_star", range(n+1), lowBound=0)
p_vars = pulp.LpVariable.dicts("p", range(n+1), lowBound=0)
w_vars = pulp.LpVariable.dicts("w", range(1, n+1), lowBound=0)

# 目标 最小化经费增量
model += delta_q, "Minimize_Delta_q"

# 约束条件
for i in range(1, n+1):
# 经费余额公式
model += p_vars[i] == (p_vars[i-1] + f1[i-1] - g_vars[i] + delta_q) * (1 + R1)

# 支出额限制
model += g_vars[i] <= p_vars[i-1] + f1[i-1] + delta_q

# 征地费用公式
model += w_vars[i] == g_vars[i] - u_vars[i-1] * (1 + R2) + u_vars[i]

# 围地数量递推关系
inv_c1 = 1.0 / c1[i-1]
model += S_star[i-1] + w_vars[i] * inv_c1 - S_prime[i-1] == S_star[i]

model += S_star[i] >= 0, f"Non-negative_S_star_{i}"
model += g_vars[i] >= 0, f"Non-negative_g_{i}"
model += u_vars[i] >= 0, f"Non-negative_u_{i}"
model += w_vars[i] >= 0, f"Non-negative_w_{i}"
model += p_vars[i] >= 0, f"Non-negative_p_{i}"

model += p_vars[0] == 0
model += p_vars[n] == 0
model += S_star[0] == 0
model += S_star[n] == 0
model += u_vars[0] == 0
model += u_vars[n] == 0

model.solve()

return delta_q.varValue, S1

# 遍历不同出矸率
e_values = np.arange(0.07, 0.11, 0.01) # 出矸率从 5% 到 15%
results = []

for e in e_values:
g0 = (d * e) / (1 - e)
v0 = g0 / rho

delta_q_opt, S1 = solve_model(v0, e)

results.append({
'e': round(e, 2),
'v0': round(v0, 2),
'S1': round(S1, 2),
'delta_q': round(delta_q_opt, 2)
})

print("{:<6} {:<12} {:<10} {:<12}".format("e (%)", "v0 (m³)", "S1 (亩)", "Δq (万元)"))
for res in results:
print(f"{res['e']*100:.0f}%\t{res['v0']}\t\t{res['S1']}\t\t{res['delta_q']}")

灵敏度稳健性分析

该论文讨论了地价涨幅 矸石容量 机械效率 电费增长 贷款利率变化 对经费的变化 通过调整相应的参数 发现地价涨幅 矸石容量 机械效率 电费增长对结果影响不大 这与论文结果一致 然而该论文并未讨论安息角与轨道仰角的灵敏度分析 并且论文中的安息角均取55°通过查询得知55°为最大安息角 论文中均取的最大安息角会导致模型高估了堆料稳定性 低估了实际占地面积 这会使模型不具推广性 仅适用于特定场景 模型得到的经费剩余偏低 我认为安息角的设置应该用保守的估计策略 例如40°左右 我们需要考虑到实际操作时产生的误差
我们通过分析安息角对经费和占地面积的影响得到了如下图表
photo

分析发现这个增长率蛮大的 因此 我认为在进行灵敏度分析的时候应该加上安息角的灵敏度分析

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
def run_model(c1, f1, S_given, R1=0.05, R2=0.05, n=5):
model = pulp.LpProblem("1999C_Coal_Gangue_Model", pulp.LpMaximize)

# 决策变量
g = pulp.LpVariable.dicts("g", range(n + 1), lowBound=0)
u = pulp.LpVariable.dicts("u", range(n + 1), lowBound=0)
w = pulp.LpVariable.dicts("w", range(n + 1), lowBound=0)
S_star = pulp.LpVariable.dicts("S_star", range(n + 1), lowBound=0)
p = pulp.LpVariable.dicts("p", range(n + 1), lowBound=0)

# 初始约束
model += p[0] == 0
model += S_star[0] == 0
model += u[0] == 0
model += u[n] == 0
model += S_star[n] == 0

# 目标函数
model += p[n]

# 约束条件
for i in range(1, n + 1):
model += w[i] == g[i] - u[i - 1] * (1 + R2) + u[i]
model += S_star[i] == S_star[i - 1] + w[i] * (1.0 / c1[i - 1]) - S_given[i - 1]
model += p[i] == (p[i - 1] + f1[i - 1] - g[i]) * (1 + R1)
model += g[i] <= p[i - 1] + f1[i - 1]

if i < n:
model += S_star[i] <= sum(S_given[k] for k in range(i, n))

# 求解
model.solve(pulp.PULP_CBC_CMD(msg=False))

# 返回最终经费余额
if pulp.LpStatus[model.status] == 'Optimal':
return p[n].varValue
else:
return None

import math

def compute_S_given(alpha_deg, base_S):
"""根据安息角调整土地需求"""
alpha_rad = math.radians(alpha_deg)
scale_factor = 1 / (math.tan(alpha_rad)**2)
return [s * scale_factor for s in base_S]

base_S = S_given[:n] # 原始占地需求
alphas = np.linspace(30, 55, 10) # 安息角从30°~55°
alpha_results = []

for alpha in alphas:
S_given_alpha = compute_S_given(alpha, base_S)
final_balance = run_model(c1=c1[:n], f1=f1[:n], S_given=S_given_alpha, n=n)
alpha_results.append({
"alpha": alpha,
"final_balance": final_balance,
"scale_factor": 1 / (math.tan(math.radians(alpha)) ** 2)
})

# 提取绘图数据
x_alphas = [r["alpha"] for r in alpha_results]
y_balances = [r["final_balance"] for r in alpha_results]
y_scales = [r["scale_factor"] for r in alpha_results]

# 安息角 vs 最终经费
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.plot(x_alphas, y_balances, marker='o', color='blue')
plt.title("安息角对最终经费的影响")
plt.xlabel("安息角 α (°)")
plt.ylabel("最终经费余额(万元)")
plt.grid()

# 安息角 vs 占地规模因子
plt.subplot(1, 2, 2)
plt.plot(x_alphas, y_scales, marker='x', color='red')
plt.title("安息角对占地规模的影响")
plt.xlabel("安息角 α (°)")
plt.ylabel("占地缩放因子")
plt.grid()

plt.tight_layout()
plt.show()

模型拓展

该论文讨论了仰角\beta的最优解 其最优解为29°
我认为可以在此基础上添加上仰角对占地面积的关系 再由占地面积推导出对经费的影响

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
# 原始参数
base_S = [10.66, 6.26, 5.25, 4.69, 4.31][:5] # 原始占地需求
c1 = [8, 8.8, 9.68, 10.648, 11.7128][:5] # 地价
f1 = [97.11, 92.33, 90.31, 88.68, 87.26][:5] # 每年收入
n = 5 # 年数

# 设置 β 的变化范围
beta_values = np.linspace(15, 43, 10) # 从 5° 到 25°,共 10 个点

results = []

for beta in beta_values:
S_given_beta = compute_S_given(beta, base_S)
final_balance = run_model(c1=c1, f1=f1, S_given=S_given_beta, n=n)

results.append({
"beta": beta,
"final_balance": final_balance,
"avg_land": sum(S_given_beta) / len(S_given_beta),
"total_land": sum(S_given_beta)
})

# 提取绘图数据
x_betas = [r["beta"] for r in results]
y_balance = [r["final_balance"] for r in results]
y_total_land = [r["total_land"] for r in results]

plt.plot(x_betas, y_total_land, marker='o', color='green')
plt.title("轨道仰角 vs 总占地面积")
plt.xlabel("轨道仰角 β (°)")
plt.ylabel("总占地面积(亩)")
plt.grid(True)

plt.tight_layout()
plt.show()

《煤矿安全规程》
第三百七十五条规定:采用轨道运输的巷道,其坡度不宜超过 15°
矿山机械设计手册》
一般矿用机车最大允许坡度为 15°~20° ;
超过 20° 的轨道需采取特殊措施

因此 我认为 这个最优解只是一个理想的最优解

  • Title: 关于1999C煤矸石堆积的研读报告
  • Author: 姜智浩
  • Created at : 2025-07-01 11:45:14
  • Updated at : 2025-07-02 11:10:42
  • Link: https://super-213.github.io/zhihaojiang.github.io/2025/07/01/20250701关于1999C煤矸石堆积的研读报告/
  • License: This work is licensed under CC BY-NC-SA 4.0.