logoChat01.ai

import os import sys import cvxpy as cp import nu...

تم الإنشاء في: ١٠ أبريل ٢٠٢٦

تمت الإجابة باستخدام GPT-5.4 Thinking بواسطة Chat01

سؤال

import os
import sys

import cvxpy as cp
import numpy as np
import pandas as pd

def run_standard_compliant_tuning(input_path, output_path):
# --- 1. 数据读取与预处理 ---
try:
df = pd.read_csv(input_path, encoding="gbk")
except:
df = pd.read_csv(input_path, encoding="utf-8")

text
# 提取原始测量数据 (根据你的列索引提取)
raw_10h = np.nan_to_num(df.iloc[:, 3].values)  # 原始10m轨向
raw_20h = np.nan_to_num(df.iloc[:, 5].values)  # 原始20m轨向
gauge_org = np.nan_to_num(df.iloc[:, 7].values)  # 原始轨距偏差
raw_10v = np.nan_to_num(df.iloc[:, 8].values)  # 原始10m高低
raw_20v = np.nan_to_num(df.iloc[:, 10].values)  # 原始20m高低
level_org = np.nan_to_num(df.iloc[:, 12].values)  # 原始水平偏差

n = len(df)
valid_10 = slice(8, -8)  # 10m弦有效区
valid_20 = slice(17, -17)  # 20m弦有效区
n_v10 = n - 16  # 10m有效点数

# --- 2. 变量定义 ---
v_l, v_r = cp.Variable(n), cp.Variable(n)  # 起道量 (Vertical)
h_l, h_r = cp.Variable(n), cp.Variable(n)  # 拨道量 (Horizontal)

# 松弛变量：用于处理“硬伤”点，防止因个别点超限导致无解
slack_l10h = cp.Variable(n_v10, nonneg=True)
slack_l10v = cp.Variable(n_v10, nonneg=True)

# --- 3. 矩阵算子构建 (弦测法物理模拟) ---
def build_matrix(size, step):
    main_diag = np.ones(size)
    side_diag = -0.5 * np.ones(size - step)
    return (
        np.diag(main_diag)
        + np.diag(side_diag, k=step)
        + np.diag(side_diag, k=-step)
    )

A10 = build_matrix(n, 8)
A20 = build_matrix(n, 17)

# --- 4. 物理模型关系式 ---
# 平面：拨道对指标的影响
new_gauge = gauge_org - h_l + h_r  # 规距 = 原规距 - 左拨 + 右拨
new_l10h = raw_10h - (A10 @ h_l)  # 10m轨向变化 (减号)
new_l20h = raw_20h - (A20 @ h_l)  # 20m轨向变化 (减号)

# 高程：起道对指标的影响
new_level = level_org - v_l + v_r  # 水平 = 原水平 - 左起 + 右起
new_l10v = raw_10v + (A10 @ v_l)  # 10m高低变化 (加号)
new_l20v = raw_20v + (A20 @ v_l)  # 20m高低变化 (加号)

# --- 5. 目标函数优化 (模拟人工经验) ---
def dead_zone(x, threshold=0.3):
    """死区函数：偏差在threshold以内不惩罚，模拟人工'保大放小'"""
    return cp.sum_squares(cp.pos(cp.abs(x) - threshold))

obj = cp.Minimize(
    # A. 平顺性优化 (权重最高，模拟视觉顺坡)
    3500.0 * cp.sum_squares(cp.diff(h_l, 2))
    + 3500.0 * cp.sum_squares(cp.diff(v_l, 2))
    # B. 10m指标优化 (使用Huber函数平衡'精准'与'顺滑')
    + 1200.0 * cp.sum(cp.huber(new_l10h[valid_10], 0.6))
    + 1200.0 * cp.sum(cp.huber(new_l10v[valid_10], 0.6))
    # C. 20m长波优化 (防止长波起伏)
    + 1000.0 * cp.sum_squares(new_l20h[valid_20])
    + 1000.0 * cp.sum_squares(new_l20v[valid_20])
    # D. 基础指标优化 (轨距、水平)
    + 600.0 * cp.sum_squares(new_gauge)
    + 600.0 * cp.sum_squares(new_level)
    # E. 惩罚项与约束模拟
    + 1000.0 * dead_zone(new_l10h[valid_10], 0.3)  # 0.3mm死区
    + 15000.0 * (cp.sum(slack_l10h) + cp.sum(slack_l10v))  # 严惩超限松弛
    + 10.0 * cp.sum_squares(h_l)  # 最小拨道量原则
    + 8.0 * cp.sum(v_l)  # 最小起道原则 (只垫不抽基准)
    # F. 边界软锁定 (给头尾一个拉向0的力，而不是死锁)
    + 800.0 * (cp.sum_squares(h_l[:17]) + cp.sum_squares(h_l[-17:]))
    + 800.0 * (cp.sum_squares(v_l[:17]) + cp.sum_squares(v_l[-17:]))
)

# --- 6. 约束条件 ---
constraints = [
    # 起拨道量程限制
    v_l >= 0,
    v_l <= 15,  # 只起不拔原则，上限15mm
    v_r >= 0,
    v_r <= 15,
    h_l >= -5,
    h_l <= 5,  # 拨道范围限制
    h_r >= -5,
    h_r <= 5,
    # 硬指标限制 (配合松弛变量)
    new_gauge >= -1.0,
    new_gauge <= 1.0,  # 轨距必须严格保1mm
    new_level >= -1.2,
    new_level <= 1.2,  # 水平保1.2mm
    new_l10h[valid_10] >= -1.2 - slack_l10h,
    new_l10h[valid_10] <= 1.2 + slack_l10h,
    new_l10v[valid_10] >= -1.2 - slack_l10v,
    new_l10v[valid_10] <= 1.2 + slack_l10v,
    # 几何安全约束
    cp.abs(new_level[5:] - new_level[:-5]) <= 1.5,  # 三角坑限值
    cp.abs(cp.diff(h_l)) <= 0.4,  # 相邻点拨道突变限制
    cp.abs(cp.diff(v_l)) <= 0.4,  # 相邻点起道突变限制
]

# --- 7. 求解运算 ---
print(f"正在进行精调优化 (当前点数: {n})...")
prob = cp.Problem(obj, constraints)
# OSQP对大规模带状矩阵鲁棒性最好
prob.solve(
    solver=cp.OSQP, verbose=False, max_iter=70000, eps_abs=1e-3, eps_rel=1e-3
)

if h_l.value is None:
    print("错误：无法找到可行解。请尝试放宽量程或检查原始数据。")
    return

# --- 8. 结果整理与导出 (0.5mm精度) ---
df.iloc[:, 13] = np.round(h_l.value * 2) / 2  # 左股拨道
df.iloc[:, 14] = np.round(h_r.value * 2) / 2  # 右股拨道

# 高程起道量 (Vertical)
df.iloc[:, 15] = np.round(v_l.value * 2) / 2  # 左股起道 (只垫不抽)
df.iloc[:, 16] = np.round(v_r.value * 2) / 2  # 右股起道 (只垫不抽)

df.to_csv(output_path, index=False, encoding="gbk")

# --- 9. 打印自检报告 ---
print("\n" + "=" * 45)
print("精调方案合规性报告:")
print(f"最大轨距偏差: {np.max(np.abs(new_gauge.value)):.2f} mm")
print(f"最大水平偏差: {np.max(np.abs(new_level.value)):.2f} mm")
print(f"最大10m轨向: {np.max(np.abs(new_l10h.value[valid_10])):.2f} mm")
print(f"最大10m高低: {np.max(np.abs(new_l10v.value[valid_10])):.2f} mm")
print(
    f"最大三角坑:   {np.max(np.abs(new_level.value[5:] - new_level.value[:-5])):.2f} mm"
)
print("-" * 45)
print(
    f"10m高低 TQI: {np.std(raw_10v[valid_10]):.3f} -> {np.std(new_l10v.value[valid_10]):.3f}"
)
print(
    f"20m高低 TQI: {np.std(raw_20v[valid_20]):.3f} -> {np.std(new_l20v.value[valid_20]):.3f}"
)
print("=" * 45)

if name == "main":
input_file = ""

text
# 1. 首先检查是否有文件拖拽进来 (命令行参数)
if len(sys.argv) > 1:
    input_file = sys.argv[1]
else:
    # 2. 如果没有拖拽文件，则进入手动输入模式
    print("--- 铁路轨道精调工具 ---")
    print("提示：你可以直接将 CSV 文件拖入此窗口，然后按回车")
    input_file = input("请输入或拖入要处理的 CSV 文件路径: ").strip()

    # 去除 Windows 路径可能带有的引号 (例如 "D:\data.csv" -> D:\data.csv)
    input_file = input_file.replace('"', "").replace("'", "")

# 3. 开始逻辑处理
if input_file:
    if os.path.exists(input_file) and input_file.lower().endswith(".csv"):
        base, ext = os.path.splitext(input_file)
        output_file = f"{base}_优化方案{ext}"

        try:
            run_standard_compliant_tuning(input_file, output_file)
            print(f"\n[成功] 方案已生成：{os.path.basename(output_file)}")
        except Exception as e:
            print(f"\n[错误] 运行崩溃，详细信息: {e}")
    else:
        print(f"\n[错误] 找不到文件或文件格式不是 CSV: {input_file}")
else:
    print("\n[错误] 未输入任何有效路径。")

print("\n" + "-" * 45)
input("按回车键退出程序...")

这是一段高铁轨道精调方案的计算代码，但是同一段轨道手动调整可以将轨向10米弦和20米弦都控制在0.5以内，这段代码只能控制在2以内，为什么

الفكرة

我先从建模本身看问题，不是先看求解器。直觉上，这段代码更像“保守顺坡 + 少动道”的优化器，而不是“把 10m/20m 弦值压到极限”的模型，所以它很容易停在 1 到 2 mm，而不是逼近 0.5 mm。

Thought for 1m 10s

الإجابة

不是求解器“算不动”，而是这段代码从建模目标上就没有在追求“10m弦、20m弦都压到 0.5 以内”。
它追求的是另一件事：少动道、顺坡、保轨距/水平、避免突变，在这个前提下尽量把指标压低。这和现场老师傅“为了把弦值打到很小，愿意更积极地分配起拨道量”的策略，不是一回事。

我直接说结论：

1. 你的 10m / 20m 指标，本质上是“软优化”，不是“硬达标”

你代码里真正的硬约束只有这些：

python
new_gauge >= -1.0
new_gauge <= 1.0
new_level >= -1.2
new_level <= 1.2

new_l10h[valid_10] >= -1.2 - slack_l10h
new_l10h[valid_10] <= 1.2 + slack_l10h
new_l10v[valid_10] >= -1.2 - slack_l10v
new_l10v[valid_10] <= 1.2 + slack_l10v

注意两点：

第一，10m 只要求到 1.2，而且还允许松弛

也就是说，模型从来没被要求“必须做到 0.5”。

第二，20m 根本没有硬约束

20m 只有这一项：

python
1000.0 * cp.sum_squares(new_l20h[valid_20])
1000.0 * cp.sum_squares(new_l20v[valid_20])

这是“尽量小”，不是“必须 ≤ 0.5”。
所以它完全可能给出一个你看起来“还行但不够狠”的结果，比如 1.2、1.5、2.0，因为模型觉得：

  • 再继续压 20m
  • 会破坏顺坡
  • 会增加起拨道量
  • 会引起相邻点突变
  • 会让头尾不自然

于是它就停了。

所以第一原因最核心：你没有把 0.5 当成约束目标。

2. 你的目标函数里，“顺坡”和“少动道”的优先级，比“把弦值压到极限”更高

看这一段：

python
3500.0 * cp.sum_squares(cp.diff(h_l, 2))
3500.0 * cp.sum_squares(cp.diff(v_l, 2))

这是二阶差分平顺性，权重非常大。
它会强烈偏好“非常顺”的起拨道曲线。

再看：

python
10.0 * cp.sum_squares(h_l)
8.0 * cp.sum(v_l)
800.0 * (cp.sum_squares(h_l[:17]) + cp.sum_squares(h_l[-17:]))
800.0 * (cp.sum_squares(v_l[:17]) + cp.sum_squares(v_l[-17:]))

这几项又在强烈表达：

  • 拨道量越小越好
  • 起道量越少越好
  • 头尾尽量别动

这就导致模型在做一种典型选择：

“我宁愿保守一点、顺一点、少动一点，也不愿为了把 10m/20m 从 1.2 再压到 0.5 去做更激进的修正。”

而人工精调恰恰常常相反：

“先把关键弦值压下去，再凭经验把动作分摊得顺一点。”

也就是说，人工是“指标优先，顺坡次之”，
而你这段代码是**“顺坡优先，指标次之”**。

3. 你给了很多“手脚镣铐”，可行域被压得很小

这些约束很硬：

python
v_l >= 0
v_l <= 15
v_r >= 0
v_r <= 15
h_l >= -5
h_l <= 5
h_r >= -5
h_r <= 5

cp.abs(new_level[5:] - new_level[:-5]) <= 1.5
cp.abs(cp.diff(h_l)) <= 0.4
cp.abs(cp.diff(v_l)) <= 0.4

尤其这两个很关键：

python
cp.abs(cp.diff(h_l)) <= 0.4
cp.abs(cp.diff(v_l)) <= 0.4

这相当于规定相邻点调整变化不能超过 0.4 mm。
如果原始波形有局部尖峰，想把 10m/20m 弦迅速压到 0.5，往往需要更灵活的局部修正。你这里把它锁得太死了。

人工为什么能更好？通常不是因为“人工比优化器聪明”，而是因为人工实际上在做这些事：

  • 某些位置敢多动一点
  • 某些位置敢让梯度陡一点
  • 左右股配合着动
  • 先处理主峰，再回头顺坡

而你的模型不太允许这种策略。

4. 物理模型很可能比现场真实调整自由度更少

你现在的弦值模型只让 左股 去影响轨向/高低：

python
new_l10h = raw_10h - (A10 @ h_l)
new_l20h = raw_20h - (A20 @ h_l)

new_l10v = raw_10v + (A10 @ v_l)
new_l20v = raw_20v + (A20 @ v_l)

右股 h_r、v_r 只参与了轨距和水平：

python
new_gauge = gauge_org - h_l + h_r
new_level = level_org - v_l + v_r

这意味着：

  • 轨向长短波全靠左股修
  • 高低长短波全靠左股修
  • 右股只是“配平轨距/水平”

如果现场人工在同一段轨道上是左右股联动来消波，那你的模型自由度就少了一半，自然不容易压得更低。

更重要的是：
如果你的原始 raw_10h/raw_20h/raw_10v/raw_20v 不是“纯左股指标”，而是检测车定义下的某种综合指标、基准线指标，
那你现在这个映射关系本身就不对，优化再努力也只能得到“假最优”。

这类问题在工程里特别常见：
不是优化不行，而是“指标—调整量”的传递模型写错了或写简化了。

5. 20m 弦的矩阵步长很可疑，可能就没按你的采样间距建对

你写的是：

python
A10 = build_matrix(n, 8)
A20 = build_matrix(n, 17)
valid_10 = slice(8, -8)
valid_20 = slice(17, -17)

这隐含一个前提：
采样间距 × 2 × step = 弦长。

比如如果点距是 0.625 m：

  • 10m 半弦对应 step = 8 没问题
  • 20m 半弦通常应接近 16，不是 17

你这里 17 很可能已经不是“20m弦”的物理算子了。

如果步长错了，会出现一个非常典型的现象：

  • 你以为你在压 20m 弦
  • 实际上你在压一个“17步离散算子”对应的东西
  • 所以和现场人工或检测系统报表的 20m 指标对不上

这个问题要重点核对：

  • 采样间距到底是多少
  • 10m/20m弦对应的半跨点数到底是多少
  • 检测系统的弦测公式是否真是
    y_i - 0.5(y_{i-k}+y_{i+k})

只要这里不一致，优化结果就不可能和人工经验完全对齐。

6. 你把 NaN 直接变成 0，会严重污染波形

python
raw_10h = np.nan_to_num(df.iloc[:, 3].values)

这类处理在几何数据里很危险。
因为 NaN 往往意味着：

  • 缺测
  • 无效区
  • 头尾不可信
  • 设备跳点

你直接补成 0，相当于在波形里硬塞了“理想点”，优化器会误以为那里本来很平顺，从而扭曲整体分配。

人工在现场不会把坏点当成 0 来算。
这也会让模型表现得比人工差。

7. 这段代码的思路更像“合规修整”，不是“极限精品打磨”

从命名其实已经暴露了：

python
run_standard_compliant_tuning

它在做的是：

  • 满足轨距/水平等安全指标
  • 尽量把 10m、20m压低
  • 同时保持非常顺
  • 同时尽量少动
  • 同时头尾别乱
  • 同时限制相邻突变

这是一个工程上偏稳健、偏保守的优化器。
而你说的人工效果“10m 和 20m 都压到 0.5 以内”，更像是：

  • 指标优先级更高
  • 容许更激进的局部修正
  • 对波形有经验性分段处理
  • 可能是分阶段、多轮微调
  • 不一定一次全局凸优化就能复现

最关键的一句话

如果人工同段确实能做到 0.5，而模型做不到，通常不是“求解器差”，而是“你的目标函数、约束和物理映射，并没有表达出人工的调法”。

你这段代码里最可能的 5 个主因，按重要性排序

第一名：20m 没有硬约束，10m 也没要求 0.5

所以模型根本没被要求做到你想要的结果。

第二名：平顺性和少动道权重过高

模型宁愿顺一点、少动一点，也不愿把弦值压到极限。

第三名：相邻点变化限制太死

abs(diff)<=0.4 很容易把局部修正能力锁死。

第四名：左右股自由度没有真正用于修正弦值

长短波几乎都压在左股上做，物理上太保守。

第五名：20m 弦离散算子可能建错

step=17 很值得怀疑。

你应该怎么改，才更接近人工效果

改法一：先把“0.5以内”写成硬目标，而不是软惩罚

比如直接加：

python
cp.abs(new_l10h[valid_10]) <= 0.5 + s10h
cp.abs(new_l10v[valid_10]) <= 0.5 + s10v
cp.abs(new_l20h[valid_20]) <= 0.5 + s20h
cp.abs(new_l20v[valid_20]) <= 0.5 + s20v

然后目标函数第一优先级先最小化这些 slack。
这样模型才会真的“拼命往 0.5 打”。

改法二：改成“两阶段优化”

这比你现在这种“一锅炖”更像人工。

第 1 阶段

先最小化：

  • 10m/20m 超 0.5 的总超限量
  • 或者最小化最大弦值

第 2 阶段

在第 1 阶段最优值固定后，再最小化：

  • 二阶差分
  • 起拨道量
  • 头尾锁定

这样就不会出现“顺坡把指标让掉”的问题。

改法三：放松这些约束试验一下

重点试这几个：

  • abs(diff(h_l)) <= 0.4 改成 0.6 或 0.8
  • abs(diff(v_l)) <= 0.4 改成 0.6 或 0.8
  • 头尾软锁定权重 800 降低
  • 平顺性权重 3500 降低

你会很快看到 10m/20m 指标显著下降。
这能直接验证“不是求解器不行，而是你给它绑太死”。

改法四：核对 20m 的离散步长

这个必须检查：

  • 点距是多少
  • 10m 对应 step 应该是多少
  • 20m 对应 step 应该是多少

如果点距不是你假设的那个，A10/A20 全都要重建。

改法五：不要把 NaN 补成 0

应该：

  • 删掉无效点
  • 或做插值
  • 或在优化时掩码掉无效区

否则波形会被你自己“造假”污染。

一个很直白的判断方法

你现在可以做个实验，不改物理模型，只改目标优先级：

  1. 把 20m 也加硬约束或强 slack
  2. 把平顺性权重从 3500 降到 500~1000
  3. 把 abs(diff)<=0.4 放宽到 0.6
  4. 去掉头尾 800 的强拉零
  5. 把 10m/20m 目标改成先最小化最大值

如果这样一改，结果立刻从 2 mm 掉到接近 0.5～0.8，说明问题就不是数据，而是策略错位。

如果这样改了还是不行，那就该重点查：

  • 弦测矩阵是否正确
  • 列索引是否取对
  • 原始指标到底对应左股、右股还是中心线
  • 手工方案是否用了你模型里没有的自由度

最后一句话

这段代码“做不到人工那样”，不是因为它弱，而是因为它被你训练成了一个：

优先保守、优先顺坡、优先少动、优先合规 的优化器，

而不是一个：

优先把 10m/20m 弦值狠狠干到 0.5 以内 的优化器。

你要是愿意，我可以下一步直接把这段代码改成更接近现场人工思路的“两阶段精调版”，把 0.5 作为一级目标写进去。

سؤال

你帮我直接把这段代码改成更接近现场人工思路的“两阶段精调版”，把 0.5 作为一级目标写进去。

الفكرة

Considering solver options and adding robustness

For solving Stage 1 and 2, I’ll stick with OSQP, as it’s suitable for quadratic programming with absolute constraints, but I’ll make sure to include fallback solvers like ECOS/SCS for robustness. I’ll also add handling for infeasible solutions and print out Stage 1 and Stage 2 results. I’ll define emergency hard caps as needed.

Thought for 2m 29s

الإجابة

可以，下面给你一个直接可替换的“两阶段精调版”。

核心变化只有三件事：

  1. 第一阶段只干一件事：把 10m/20m 轨向、高低 超过 0.5 mm 的部分尽可能压掉。
  2. 第二阶段再做顺坡、少动道、头尾过渡，但不允许明显破坏第一阶段成果。
  3. 把原来直接 nan_to_num(...)=0 改成了插值补点，避免坏点把波形拉假平。

直接替换版代码

python
import os
import sys

import cvxpy as cp
import numpy as np
import pandas as pd


# =========================
# 可调参数
# =========================
TARGET_LIMIT = 0.5   # 一级目标：10m / 20m 尽量压到 0.5 mm 以内
HARD_LIMIT = 1.2     # 二级硬上限：若 0.5 无法完全实现，仍尽量保在 1.2 以内
GAUGE_LIMIT = 1.0    # 轨距偏差硬限
LEVEL_LIMIT = 1.2    # 水平偏差硬限
TRIANGLE_LIMIT = 1.5 # 三角坑限值

# 采样步长：必须与你的轨检/测量点距一致
STEP_10 = 8
STEP_20 = 17  # 注意：这里一定要按你的检测系统定义核对，很多系统 20m 半弦未必是 17

MAX_LIFT = 15.0
MAX_LINE = 5.0

# 比原版稍微放开一点，更接近“先压峰再顺坡”
MAX_DIFF_H = 0.6
MAX_DIFF_V = 0.6

ROUND_TO = 0.5  # 输出量取 0.5mm 精度


def read_csv_auto(path):
    try:
        return pd.read_csv(path, encoding="gbk")
    except Exception:
        return pd.read_csv(path, encoding="utf-8")


def clean_numeric_col(df, idx):
    """读取数值列，优先做插值，不直接把 NaN 当 0。"""
    s = pd.to_numeric(df.iloc[:, idx], errors="coerce").astype(float)
    if s.notna().sum() == 0:
        raise ValueError(f"第 {idx} 列没有可用数值，无法计算。")
    s = s.interpolate(limit_direction="both")
    s = s.fillna(method="ffill").fillna(method="bfill").fillna(0.0)
    return s.to_numpy(dtype=float)


def build_matrix(size, step):
    """弦测法离散算子：y_i - 0.5*(y_{i-step}+y_{i+step})"""
    main_diag = np.ones(size)
    side_diag = -0.5 * np.ones(size - step)
    return (
        np.diag(main_diag)
        + np.diag(side_diag, k=step)
        + np.diag(side_diag, k=-step)
    )


def round_half(x):
    return np.round(x / ROUND_TO) * ROUND_TO


def solve_with_fallback(prob, stage_name):
    """优先 OSQP，失败则回退到 SCS。"""
    solvers = [
        ("OSQP", cp.OSQP, dict(verbose=False, warm_start=True, max_iter=80000, eps_abs=1e-4, eps_rel=1e-4)),
        ("SCS", cp.SCS, dict(verbose=False, warm_start=True, max_iters=30000, eps=1e-4)),
    ]

    last_err = None
    for solver_name, solver, kwargs in solvers:
        try:
            prob.solve(solver=solver, **kwargs)
            if prob.status in ("optimal", "optimal_inaccurate"):
                print(f"{stage_name} 求解成功，solver={solver_name}, status={prob.status}")
                return
        except Exception as e:
            last_err = e

    raise RuntimeError(f"{stage_name} 求解失败。最后错误: {last_err}")


def run_two_stage_manual_style_tuning(input_path, output_path):
    # =========================
    # 1. 读数
    # =========================
    df = read_csv_auto(input_path)

    raw_10h = clean_numeric_col(df, 3)   # 原始10m轨向
    raw_20h = clean_numeric_col(df, 5)   # 原始20m轨向
    gauge_org = clean_numeric_col(df, 7) # 原始轨距偏差

    raw_10v = clean_numeric_col(df, 8)   # 原始10m高低
    raw_20v = clean_numeric_col(df, 10)  # 原始20m高低
    level_org = clean_numeric_col(df, 12)# 原始水平偏差

    n = len(df)
    if n <= 2 * STEP_20:
        raise ValueError(f"数据点数不足：当前 {n} 点，不足以计算 20m 弦有效区。")

    valid_10 = slice(STEP_10, -STEP_10)
    valid_20 = slice(STEP_20, -STEP_20)

    m10 = n - 2 * STEP_10
    m20 = n - 2 * STEP_20

    A10 = build_matrix(n, STEP_10)
    A20 = build_matrix(n, STEP_20)

    # =========================
    # 2. 变量
    # =========================
    v_l = cp.Variable(n)  # 左股起道
    v_r = cp.Variable(n)  # 右股起道
    h_l = cp.Variable(n)  # 左股拨道
    h_r = cp.Variable(n)  # 右股拨道

    # =========================
    # 3. 物理模型
    # =========================
    # 平面
    new_gauge = gauge_org - h_l + h_r
    new_l10h = raw_10h - (A10 @ h_l)
    new_l20h = raw_20h - (A20 @ h_l)

    # 高程
    new_level = level_org - v_l + v_r
    new_l10v = raw_10v + (A10 @ v_l)
    new_l20v = raw_20v + (A20 @ v_l)

    # =========================
    # 4. 基础硬约束
    # =========================
    base_constraints = [
        # 起拨道量程
        v_l >= 0,
        v_l <= MAX_LIFT,
        v_r >= 0,
        v_r <= MAX_LIFT,

        h_l >= -MAX_LINE,
        h_l <= MAX_LINE,
        h_r >= -MAX_LINE,
        h_r <= MAX_LINE,

        # 安全硬指标
        new_gauge >= -GAUGE_LIMIT,
        new_gauge <= GAUGE_LIMIT,
        new_level >= -LEVEL_LIMIT,
        new_level <= LEVEL_LIMIT,

        # 几何约束
        cp.abs(new_level[5:] - new_level[:-5]) <= TRIANGLE_LIMIT,

        # 比原版略放松，更接近现场“先压峰”
        cp.abs(cp.diff(h_l)) <= MAX_DIFF_H,
        cp.abs(cp.diff(v_l)) <= MAX_DIFF_V,
        cp.abs(cp.diff(h_r)) <= MAX_DIFF_H,
        cp.abs(cp.diff(v_r)) <= MAX_DIFF_V,
    ]

    # =========================
    # 5. 第一阶段：先压 0.5 目标
    # =========================
    # exceed_* 表示“超出 0.5 的部分”
    exceed_10h = cp.Variable(m10, nonneg=True)
    exceed_10v = cp.Variable(m10, nonneg=True)
    exceed_20h = cp.Variable(m20, nonneg=True)
    exceed_20v = cp.Variable(m20, nonneg=True)

    # emergency_* 表示“超出 1.2 的部分”，仅作为兜底，防止局部异常失控
    emergency_10h = cp.Variable(m10, nonneg=True)
    emergency_10v = cp.Variable(m10, nonneg=True)
    emergency_20h = cp.Variable(m20, nonneg=True)
    emergency_20v = cp.Variable(m20, nonneg=True)

    peak_exceed = cp.Variable(nonneg=True)  # 最大超 0.5 量（峰值）
    total_exceed = (
        cp.sum(exceed_10h) + cp.sum(exceed_10v)
        + cp.sum(exceed_20h) + cp.sum(exceed_20v)
    )
    total_emergency = (
        cp.sum(emergency_10h) + cp.sum(emergency_10v)
        + cp.sum(emergency_20h) + cp.sum(emergency_20v)
    )

    stage1_constraints = base_constraints + [
        # 超 0.5 的定义：exceed >= |x| - 0.5
        exceed_10h >= cp.abs(new_l10h[valid_10]) - TARGET_LIMIT,
        exceed_10v >= cp.abs(new_l10v[valid_10]) - TARGET_LIMIT,
        exceed_20h >= cp.abs(new_l20h[valid_20]) - TARGET_LIMIT,
        exceed_20v >= cp.abs(new_l20v[valid_20]) - TARGET_LIMIT,

        # 最大超限峰值
        exceed_10h <= peak_exceed,
        exceed_10v <= peak_exceed,
        exceed_20h <= peak_exceed,
        exceed_20v <= peak_exceed,

        # 兜底硬上限：尽量不超过 1.2，超了就付极高代价
        cp.abs(new_l10h[valid_10]) <= HARD_LIMIT + emergency_10h,
        cp.abs(new_l10v[valid_10]) <= HARD_LIMIT + emergency_10v,
        cp.abs(new_l20h[valid_20]) <= HARD_LIMIT + emergency_20h,
        cp.abs(new_l20v[valid_20]) <= HARD_LIMIT + emergency_20v,
    ]

    # 第一阶段：压峰优先，总超限次之，兜底超限严惩
    stage1_obj = cp.Minimize(
        1_000_000.0 * peak_exceed
        + 20_000.0 * (cp.sum(exceed_10h) + cp.sum(exceed_10v))
        + 12_000.0 * (cp.sum(exceed_20h) + cp.sum(exceed_20v))
        + 300_000.0 * total_emergency
        # 很小的正则，防止第一阶段解漂移
        + 20.0 * cp.sum_squares(h_l)
        + 20.0 * cp.sum_squares(h_r)
        + 8.0 * (cp.sum(v_l) + cp.sum(v_r))
        + 50.0 * cp.sum_squares(cp.diff(h_l, 2))
        + 50.0 * cp.sum_squares(cp.diff(v_l, 2))
    )

    print(f"正在执行第一阶段（先压 0.5 目标），点数: {n}")
    prob1 = cp.Problem(stage1_obj, stage1_constraints)
    solve_with_fallback(prob1, "第一阶段")

    if h_l.value is None:
        raise RuntimeError("第一阶段未得到可行解。")

    best_peak_exceed = float(peak_exceed.value)
    best_total_exceed = float(total_exceed.value)
    best_total_emergency = float(total_emergency.value)

    print("-" * 50)
    print("第一阶段结果（压峰结果）:")
    print(f"最大超 0.5 量: {best_peak_exceed:.4f} mm")
    print(f"总超 0.5 量:   {best_total_exceed:.4f}")
    print(f"超 1.2 总量:   {best_total_emergency:.6f}")
    print("-" * 50)

    # =========================
    # 6. 第二阶段：在不破坏第一阶段成果下，做顺坡/少动/过渡
    # =========================
    lock_tol_peak = 1e-3
    lock_tol_sum = 1e-2
    lock_tol_emg = 1e-6

    stage2_constraints = stage1_constraints + [
        peak_exceed <= best_peak_exceed + lock_tol_peak,
        total_exceed <= best_total_exceed + lock_tol_sum,
        total_emergency <= best_total_emergency + lock_tol_emg,
    ]

    stage2_obj = cp.Minimize(
        # A. 平顺性：主目标
        3200.0 * cp.sum_squares(cp.diff(h_l, 2))
        + 3200.0 * cp.sum_squares(cp.diff(v_l, 2))
        + 1200.0 * cp.sum_squares(cp.diff(h_r, 2))
        + 1200.0 * cp.sum_squares(cp.diff(v_r, 2))

        # B. 让已压下来的 10m/20m 进一步更圆顺
        + 260.0 * cp.sum_squares(new_l10h[valid_10])
        + 260.0 * cp.sum_squares(new_l10v[valid_10])
        + 180.0 * cp.sum_squares(new_l20h[valid_20])
        + 180.0 * cp.sum_squares(new_l20v[valid_20])

        # C. 轨距/水平同步优化
        + 700.0 * cp.sum_squares(new_gauge)
        + 700.0 * cp.sum_squares(new_level)

        # D. 少动道原则
        + 18.0 * cp.sum_squares(h_l)
        + 18.0 * cp.sum_squares(h_r)
        + 6.0 * (cp.sum(v_l) + cp.sum(v_r))

        # E. 头尾过渡，不像原版那样锁太死
        + 180.0 * (
            cp.sum_squares(h_l[:STEP_20]) + cp.sum_squares(h_l[-STEP_20:])
            + cp.sum_squares(v_l[:STEP_20]) + cp.sum_squares(v_l[-STEP_20:])
            + cp.sum_squares(h_r[:STEP_20]) + cp.sum_squares(h_r[-STEP_20:])
            + cp.sum_squares(v_r[:STEP_20]) + cp.sum_squares(v_r[-STEP_20:])
        )
    )

    print("正在执行第二阶段（顺坡 + 少动道 + 头尾过渡）...")
    prob2 = cp.Problem(stage2_obj, stage2_constraints)
    solve_with_fallback(prob2, "第二阶段")

    if h_l.value is None:
        raise RuntimeError("第二阶段未得到可行解。")

    # =========================
    # 7. 写回结果
    # =========================
    df.iloc[:, 13] = round_half(h_l.value)  # 左股拨道
    df.iloc[:, 14] = round_half(h_r.value)  # 右股拨道
    df.iloc[:, 15] = round_half(v_l.value)  # 左股起道
    df.iloc[:, 16] = round_half(v_r.value)  # 右股起道

    df.to_csv(output_path, index=False, encoding="gbk")

    # =========================
    # 8. 自检报告
    # =========================
    final_10h = np.abs(new_l10h.value[valid_10])
    final_20h = np.abs(new_l20h.value[valid_20])
    final_10v = np.abs(new_l10v.value[valid_10])
    final_20v = np.abs(new_l20v.value[valid_20])

    def ratio_leq(arr, threshold):
        return 100.0 * np.mean(arr <= threshold)

    print("\n" + "=" * 58)
    print("两阶段精调方案报告:")
    print(f"最大轨距偏差: {np.max(np.abs(new_gauge.value)):.3f} mm")
    print(f"最大水平偏差: {np.max(np.abs(new_level.value)):.3f} mm")
    print(f"最大10m轨向:   {np.max(final_10h):.3f} mm")
    print(f"最大20m轨向:   {np.max(final_20h):.3f} mm")
    print(f"最大10m高低:   {np.max(final_10v):.3f} mm")
    print(f"最大20m高低:   {np.max(final_20v):.3f} mm")
    print(f"最大三角坑:    {np.max(np.abs(new_level.value[5:] - new_level.value[:-5])):.3f} mm")
    print("-" * 58)
    print(f"10m轨向 ≤0.5 占比: {ratio_leq(final_10h, 0.5):6.2f}%")
    print(f"20m轨向 ≤0.5 占比: {ratio_leq(final_20h, 0.5):6.2f}%")
    print(f"10m高低 ≤0.5 占比: {ratio_leq(final_10v, 0.5):6.2f}%")
    print(f"20m高低 ≤0.5 占比: {ratio_leq(final_20v, 0.5):6.2f}%")
    print("-" * 58)
    print(
        f"10m高低 TQI: {np.std(raw_10v[valid_10]):.3f} -> {np.std(new_l10v.value[valid_10]):.3f}"
    )
    print(
        f"20m高低 TQI: {np.std(raw_20v[valid_20]):.3f} -> {np.std(new_l20v.value[valid_20]):.3f}"
    )
    print(
        f"10m轨向 TQI: {np.std(raw_10h[valid_10]):.3f} -> {np.std(new_l10h.value[valid_10]):.3f}"
    )
    print(
        f"20m轨向 TQI: {np.std(raw_20h[valid_20]):.3f} -> {np.std(new_l20h.value[valid_20]):.3f}"
    )
    print("=" * 58)


if __name__ == "__main__":
    input_file = ""

    if len(sys.argv) > 1:
        input_file = sys.argv[1]
    else:
        print("--- 铁路轨道两阶段精调工具 ---")
        print("提示：你可以直接将 CSV 文件拖入此窗口，然后按回车")
        input_file = input("请输入或拖入要处理的 CSV 文件路径: ").strip()
        input_file = input_file.replace('"', "").replace("'", "")

    if input_file:
        if os.path.exists(input_file) and input_file.lower().endswith(".csv"):
            base, ext = os.path.splitext(input_file)
            output_file = f"{base}_两阶段优化方案{ext}"

            try:
                run_two_stage_manual_style_tuning(input_file, output_file)
                print(f"\n[成功] 方案已生成：{os.path.basename(output_file)}")
            except Exception as e:
                print(f"\n[错误] 运行失败，详细信息: {e}")
        else:
            print(f"\n[错误] 找不到文件或文件格式不是 CSV: {input_file}")
    else:
        print("\n[错误] 未输入任何有效路径。")

    print("\n" + "-" * 45)
    input("按回车键退出程序...")

这版和你原版最关键的差别

原版是：

  • 顺坡优先
  • 少动道优先
  • 10m/20m 只是“尽量优化”

这版改成：

  • 第一阶段先看 0.5
  • 第二阶段才看顺坡和少动道

这就更接近现场人工那种“先打主峰，再做圆顺”的路子。

你拿去跑之前，最应该先检查的两个点

1. STEP_20 = 17 到底对不对

这个非常关键。
如果你的采样间距不是恰好对应 17 个点半弦，那 20m 模型就会偏。
这会直接导致“人工能压到 0.5，程序压不下去”。

2. 10m/20m 指标是不是确实只对应左股

你现在这个物理模型还是沿用了你原来的假设：

  • 轨向/高低长短波主要由左股变化承担
  • 右股主要负责轨距/水平配平

如果你的现场检测指标其实是中心线、或左右股综合定义，那还要再改模型。

这版如果还压不到 0.5，下一步该怎么调

优先动这几个参数：

  • MAX_DIFF_H = 0.6 可试到 0.8
  • MAX_DIFF_V = 0.6 可试到 0.8
  • 第二阶段里平顺性权重 3200 可降到 2200
  • 头尾过渡权重 180 可降到 80

如果你愿意，我下一步可以继续帮你做一版：

“分段主峰优先 + 左右股联动分摊 + 自动参数扫描”的现场强化版。

سؤال

直接把这段代码改成更接近现场人工思路的“两阶段精调版”，把 0.5 作为一级目标写进去。

الفكرة

Deciding on the code response format

The user asked for code modification again, and it seems like they just want a final, concise version with the code block. I'll focus on this direct request and provide the modified code. No need for extra container or execution, just a quick response with key changes included.

Providing the revised code in Chinese

The user likely wants the revised code, so I should provide the full, updated script, ensuring it's polished and better than the initial version. I'll enhance it with lexicographic optimization and mention it's directly replaceable. I'll keep original structure and comments, focusing on cvxpy compatibility and correct expressions.

Thought for 19s

الإجابة

下面这版可以直接替换。思路就是你要的：

  • 第一阶段：先把 10m/20m 轨向、高低 超过 0.5 mm 的部分尽可能压掉
  • 第二阶段：在不明显破坏第一阶段成果的前提下，再做顺坡、少动道、头尾过渡

我尽量保留了你原来的列位和使用方式。

python
import os
import sys

import cvxpy as cp
import numpy as np
import pandas as pd


# =========================
# 参数区：按现场情况微调
# =========================
TARGET_LIMIT = 0.5     # 一级目标：10m/20m 尽量压到 0.5 以内
BACKUP_LIMIT = 1.2     # 兜底硬限
GAUGE_LIMIT = 1.0      # 轨距偏差硬限
LEVEL_LIMIT = 1.2      # 水平偏差硬限
TRIANGLE_LIMIT = 1.5   # 三角坑限值

STEP_10 = 8
STEP_20 = 17  # 这个必须和你的点距/检测定义核对

MAX_LIFT = 15.0        # 起道上限
MAX_LINE = 5.0         # 拨道上限

MAX_DIFF_H = 0.6       # 相邻点拨道变化限值（比原版略放松）
MAX_DIFF_V = 0.6       # 相邻点起道变化限值（比原版略放松）

ROUND_PRECISION = 0.5  # 输出按 0.5mm 取整


def read_csv_auto(path):
    try:
        return pd.read_csv(path, encoding="gbk")
    except Exception:
        return pd.read_csv(path, encoding="utf-8")


def get_numeric_column(df, idx):
    """
    数值列清洗：
    1. 非数值转 NaN
    2. 用插值补缺，不直接把 NaN 当 0
    """
    s = pd.to_numeric(df.iloc[:, idx], errors="coerce").astype(float)
    if s.notna().sum() == 0:
        raise ValueError(f"第 {idx} 列没有有效数值。")
    s = s.interpolate(limit_direction="both")
    s = s.ffill().bfill().fillna(0.0)
    return s.to_numpy(dtype=float)


def build_matrix(size, step):
    """
    弦测法离散算子：
    y_i - 0.5 * (y_{i-step} + y_{i+step})
    """
    main_diag = np.ones(size)
    side_diag = -0.5 * np.ones(size - step)
    return (
        np.diag(main_diag)
        + np.diag(side_diag, k=step)
        + np.diag(side_diag, k=-step)
    )


def round_to_step(x, step=0.5):
    return np.round(x / step) * step


def solve_problem(prob, stage_name):
    """
    优先 OSQP，失败时回退 SCS
    """
    tried = []

    try:
        prob.solve(
            solver=cp.OSQP,
            verbose=False,
            warm_start=True,
            max_iter=80000,
            eps_abs=1e-4,
            eps_rel=1e-4,
        )
        tried.append(("OSQP", prob.status))
        if prob.status in ("optimal", "optimal_inaccurate"):
            print(f"{stage_name} 求解成功: OSQP / {prob.status}")
            return
    except Exception as e:
        tried.append(("OSQP", str(e)))

    try:
        prob.solve(
            solver=cp.SCS,
            verbose=False,
            warm_start=True,
            max_iters=30000,
            eps=1e-4,
        )
        tried.append(("SCS", prob.status))
        if prob.status in ("optimal", "optimal_inaccurate"):
            print(f"{stage_name} 求解成功: SCS / {prob.status}")
            return
    except Exception as e:
        tried.append(("SCS", str(e)))

    raise RuntimeError(f"{stage_name} 求解失败，尝试记录: {tried}")


def run_two_stage_manual_style_tuning(input_path, output_path):
    # =========================
    # 1. 数据读取
    # =========================
    df = read_csv_auto(input_path)

    raw_10h = get_numeric_column(df, 3)    # 原始10m轨向
    raw_20h = get_numeric_column(df, 5)    # 原始20m轨向
    gauge_org = get_numeric_column(df, 7)  # 原始轨距偏差

    raw_10v = get_numeric_column(df, 8)    # 原始10m高低
    raw_20v = get_numeric_column(df, 10)   # 原始20m高低
    level_org = get_numeric_column(df, 12) # 原始水平偏差

    n = len(df)
    if n <= 2 * STEP_20:
        raise ValueError(f"数据点数不足，当前 {n} 点，无法形成 20m 有效区。")

    valid_10 = slice(STEP_10, -STEP_10)
    valid_20 = slice(STEP_20, -STEP_20)

    m10 = n - 2 * STEP_10
    m20 = n - 2 * STEP_20

    A10 = build_matrix(n, STEP_10)
    A20 = build_matrix(n, STEP_20)

    # =========================
    # 2. 决策变量
    # =========================
    v_l = cp.Variable(n)  # 左股起道
    v_r = cp.Variable(n)  # 右股起道
    h_l = cp.Variable(n)  # 左股拨道
    h_r = cp.Variable(n)  # 右股拨道

    # =========================
    # 3. 物理模型
    # =========================
    # 平面
    new_gauge = gauge_org - h_l + h_r
    new_l10h = raw_10h - (A10 @ h_l)
    new_l20h = raw_20h - (A20 @ h_l)

    # 高程
    new_level = level_org - v_l + v_r
    new_l10v = raw_10v + (A10 @ v_l)
    new_l20v = raw_20v + (A20 @ v_l)

    # =========================
    # 4. 公共硬约束
    # =========================
    base_constraints = [
        # 量程约束
        v_l >= 0,
        v_l <= MAX_LIFT,
        v_r >= 0,
        v_r <= MAX_LIFT,
        h_l >= -MAX_LINE,
        h_l <= MAX_LINE,
        h_r >= -MAX_LINE,
        h_r <= MAX_LINE,

        # 安全硬约束
        new_gauge >= -GAUGE_LIMIT,
        new_gauge <= GAUGE_LIMIT,
        new_level >= -LEVEL_LIMIT,
        new_level <= LEVEL_LIMIT,

        # 三角坑
        cp.abs(new_level[5:] - new_level[:-5]) <= TRIANGLE_LIMIT,

        # 相邻点变化限制：比原版略放开，更像先压峰再顺坡
        cp.abs(cp.diff(h_l)) <= MAX_DIFF_H,
        cp.abs(cp.diff(h_r)) <= MAX_DIFF_H,
        cp.abs(cp.diff(v_l)) <= MAX_DIFF_V,
        cp.abs(cp.diff(v_r)) <= MAX_DIFF_V,
    ]

    # =========================
    # 5. 第一阶段：把 0.5 作为一级目标
    # =========================
    # exceed_* 表示 “超过 0.5 的那部分”
    exceed_10h = cp.Variable(m10, nonneg=True)
    exceed_10v = cp.Variable(m10, nonneg=True)
    exceed_20h = cp.Variable(m20, nonneg=True)
    exceed_20v = cp.Variable(m20, nonneg=True)

    # backup_* 表示 “超过 1.2 的那部分”，极重惩罚，兜底
    backup_10h = cp.Variable(m10, nonneg=True)
    backup_10v = cp.Variable(m10, nonneg=True)
    backup_20h = cp.Variable(m20, nonneg=True)
    backup_20v = cp.Variable(m20, nonneg=True)

    peak_exceed = cp.Variable(nonneg=True)

    total_exceed = (
        cp.sum(exceed_10h) + cp.sum(exceed_10v)
        + cp.sum(exceed_20h) + cp.sum(exceed_20v)
    )
    total_backup = (
        cp.sum(backup_10h) + cp.sum(backup_10v)
        + cp.sum(backup_20h) + cp.sum(backup_20v)
    )

    stage1_constraints = base_constraints + [
        # exceed >= |x| - 0.5
        exceed_10h >= cp.abs(new_l10h[valid_10]) - TARGET_LIMIT,
        exceed_10v >= cp.abs(new_l10v[valid_10]) - TARGET_LIMIT,
        exceed_20h >= cp.abs(new_l20h[valid_20]) - TARGET_LIMIT,
        exceed_20v >= cp.abs(new_l20v[valid_20]) - TARGET_LIMIT,

        # 把峰值也拉下来：相当于优先打主峰
        exceed_10h <= peak_exceed,
        exceed_10v <= peak_exceed,
        exceed_20h <= peak_exceed,
        exceed_20v <= peak_exceed,

        # 兜底：尽量不要超过 1.2
        cp.abs(new_l10h[valid_10]) <= BACKUP_LIMIT + backup_10h,
        cp.abs(new_l10v[valid_10]) <= BACKUP_LIMIT + backup_10v,
        cp.abs(new_l20h[valid_20]) <= BACKUP_LIMIT + backup_20h,
        cp.abs(new_l20v[valid_20]) <= BACKUP_LIMIT + backup_20v,
    ]

    # 第一阶段目标：
    # 1) 先尽量压低最大超 0.5 的峰
    # 2) 再压总超限量
    # 3) 超过 1.2 的部分极重惩罚
    # 4) 只给一点轻微正则，防止解飘
    stage1_obj = cp.Minimize(
        1_000_000.0 * peak_exceed
        + 20_000.0 * (cp.sum(exceed_10h) + cp.sum(exceed_10v))
        + 12_000.0 * (cp.sum(exceed_20h) + cp.sum(exceed_20v))
        + 300_000.0 * total_backup
        + 15.0 * cp.sum_squares(h_l)
        + 15.0 * cp.sum_squares(h_r)
        + 6.0 * (cp.sum(v_l) + cp.sum(v_r))
        + 30.0 * cp.sum_squares(cp.diff(h_l, 2))
        + 30.0 * cp.sum_squares(cp.diff(v_l, 2))
    )

    print(f"正在执行第一阶段：先压 10m/20m 到 0.5 附近（点数 {n}）...")
    prob1 = cp.Problem(stage1_obj, stage1_constraints)
    solve_problem(prob1, "第一阶段")

    if h_l.value is None:
        raise RuntimeError("第一阶段未找到可行解。")

    best_peak_exceed = float(peak_exceed.value)
    best_total_exceed = float(total_exceed.value)
    best_total_backup = float(total_backup.value)

    print("-" * 50)
    print("第一阶段完成：")
    print(f"最大超 0.5 量: {best_peak_exceed:.4f} mm")
    print(f"总超 0.5 量:   {best_total_exceed:.4f}")
    print(f"总超 1.2 量:   {best_total_backup:.6f}")
    print("-" * 50)

    # =========================
    # 6. 第二阶段：在不破坏第一阶段成果下，做顺坡和少动道
    # =========================
    # 给一点极小容差，避免数值误差导致不可行
    stage2_constraints = stage1_constraints + [
        peak_exceed <= best_peak_exceed + 1e-3,
        total_exceed <= best_total_exceed + 1e-2,
        total_backup <= best_total_backup + 1e-6,
    ]

    stage2_obj = cp.Minimize(
        # A. 顺坡：主目标
        3200.0 * cp.sum_squares(cp.diff(h_l, 2))
        + 3200.0 * cp.sum_squares(cp.diff(v_l, 2))
        + 1000.0 * cp.sum_squares(cp.diff(h_r, 2))
        + 1000.0 * cp.sum_squares(cp.diff(v_r, 2))

        # B. 保持 10m/20m 圆顺
        + 240.0 * cp.sum_squares(new_l10h[valid_10])
        + 240.0 * cp.sum_squares(new_l10v[valid_10])
        + 160.0 * cp.sum_squares(new_l20h[valid_20])
        + 160.0 * cp.sum_squares(new_l20v[valid_20])

        # C. 轨距、水平继续优化
        + 700.0 * cp.sum_squares(new_gauge)
        + 700.0 * cp.sum_squares(new_level)

        # D. 少动道
        + 16.0 * cp.sum_squares(h_l)
        + 16.0 * cp.sum_squares(h_r)
        + 5.0 * (cp.sum(v_l) + cp.sum(v_r))

        # E. 头尾过渡，不像原版那样锁太死
        + 180.0 * (
            cp.sum_squares(h_l[:STEP_20]) + cp.sum_squares(h_l[-STEP_20:])
            + cp.sum_squares(h_r[:STEP_20]) + cp.sum_squares(h_r[-STEP_20:])
            + cp.sum_squares(v_l[:STEP_20]) + cp.sum_squares(v_l[-STEP_20:])
            + cp.sum_squares(v_r[:STEP_20]) + cp.sum_squares(v_r[-STEP_20:])
        )
    )

    print("正在执行第二阶段：顺坡、少动道、头尾过渡...")
    prob2 = cp.Problem(stage2_obj, stage2_constraints)
    solve_problem(prob2, "第二阶段")

    if h_l.value is None:
        raise RuntimeError("第二阶段未找到可行解。")

    # =========================
    # 7. 写回结果
    # =========================
    df.iloc[:, 13] = round_to_step(h_l.value, ROUND_PRECISION)  # 左股拨道
    df.iloc[:, 14] = round_to_step(h_r.value, ROUND_PRECISION)  # 右股拨道
    df.iloc[:, 15] = round_to_step(v_l.value, ROUND_PRECISION)  # 左股起道
    df.iloc[:, 16] = round_to_step(v_r.value, ROUND_PRECISION)  # 右股起道

    df.to_csv(output_path, index=False, encoding="gbk")

    # =========================
    # 8. 报告
    # =========================
    final_10h = np.abs(new_l10h.value[valid_10])
    final_20h = np.abs(new_l20h.value[valid_20])
    final_10v = np.abs(new_l10v.value[valid_10])
    final_20v = np.abs(new_l20v.value[valid_20])

    def ratio_leq(arr, th):
        return 100.0 * np.mean(arr <= th)

    print("\n" + "=" * 58)
    print("两阶段精调方案报告")
    print(f"最大轨距偏差: {np.max(np.abs(new_gauge.value)):.3f} mm")
    print(f"最大水平偏差: {np.max(np.abs(new_level.value)):.3f} mm")
    print(f"最大10m轨向:   {np.max(final_10h):.3f} mm")
    print(f"最大20m轨向:   {np.max(final_20h):.3f} mm")
    print(f"最大10m高低:   {np.max(final_10v):.3f} mm")
    print(f"最大20m高低:   {np.max(final_20v):.3f} mm")
    print(f"最大三角坑:    {np.max(np.abs(new_level.value[5:] - new_level.value[:-5])):.3f} mm")
    print("-" * 58)
    print(f"10m轨向 <= 0.5 占比: {ratio_leq(final_10h, 0.5):6.2f}%")
    print(f"20m轨向 <= 0.5 占比: {ratio_leq(final_20h, 0.5):6.2f}%")
    print(f"10m高低 <= 0.5 占比: {ratio_leq(final_10v, 0.5):6.2f}%")
    print(f"20m高低 <= 0.5 占比: {ratio_leq(final_20v, 0.5):6.2f}%")
    print("-" * 58)
    print(f"10m轨向TQI: {np.std(raw_10h[valid_10]):.3f} -> {np.std(new_l10h.value[valid_10]):.3f}")
    print(f"20m轨向TQI: {np.std(raw_20h[valid_20]):.3f} -> {np.std(new_l20h.value[valid_20]):.3f}")
    print(f"10m高低TQI: {np.std(raw_10v[valid_10]):.3f} -> {np.std(new_l10v.value[valid_10]):.3f}")
    print(f"20m高低TQI: {np.std(raw_20v[valid_20]):.3f} -> {np.std(new_l20v.value[valid_20]):.3f}")
    print("=" * 58)


if __name__ == "__main__":
    input_file = ""

    if len(sys.argv) > 1:
        input_file = sys.argv[1]
    else:
        print("--- 铁路轨道两阶段精调工具 ---")
        print("提示：你可以直接将 CSV 文件拖入此窗口，然后按回车")
        input_file = input("请输入或拖入要处理的 CSV 文件路径: ").strip()
        input_file = input_file.replace('"', "").replace("'", "")

    if input_file:
        if os.path.exists(input_file) and input_file.lower().endswith(".csv"):
            base, ext = os.path.splitext(input_file)
            output_file = f"{base}_两阶段优化方案{ext}"

            try:
                run_two_stage_manual_style_tuning(input_file, output_file)
                print(f"\n[成功] 方案已生成：{os.path.basename(output_file)}")
            except Exception as e:
                print(f"\n[错误] 运行失败，详细信息: {e}")
        else:
            print(f"\n[错误] 找不到文件或文件格式不是 CSV: {input_file}")
    else:
        print("\n[错误] 未输入任何有效路径。")

    print("\n" + "-" * 45)
    input("按回车键退出程序...")

这版已经把“0.5 作为一级目标”写进去了。实际含义是：

  • 不是简单地“尽量小”
  • 而是先最小化 超过 0.5 的那部分
  • 再在第二阶段做顺坡和少动

你跑这版时，优先核对一件事：STEP_20 = 17 是否真的对应你的 20m 弦定义。这个如果错了，20m 怎么调都容易和现场不一致。

سؤال

你调整的这版本代码运行结果比我最初发给你的那一版还要差，你有什么头绪吗

الفكرة

我先不继续“拍脑袋调权重”了，先把问题缩到几处最可能的结构性原因：弦测模型、左右股作用方式、以及第二阶段是不是把第一阶段成果又冲掉了。先查一下有没有现成的代码分析类 skill 可借用，再给你一个更靠谱的排查结论。

Thought for 10s

الفكرة

Analyzing code constraints and optimization

It looks like stage1 constraints are controlling exceed values due to peak_exceed limits. In stage2, while the exceed values aren't directly minimized since the objective's focus is on totals and bounds, these variables can increase freely as long as they stay within limits. However, slacks in the constraints might inadvertently allow larger values in stage2.

Reviewing Stage2 constraints and optimization

In Stage2, the constraints primarily lock the maximum exceed and total exceed values relative to the best peak and total from Stage1. Stage2 can adjust exceed allocations to absorb worsened values, which might lead to individual points degrading while keeping the total within bounds. This might cause further issues if Stage1’s solution wasn’t ideal.

Thought for 1m 12s

الإجابة

有，而且我现在最怀疑的不是“参数没调好”，而是我那版改法里有几处结构性问题。按优先级说：

1. 第二阶段很可能把第一阶段的成果“冲掉了”

这是我最先怀疑的。

我那版第二阶段虽然加了这些锁：

python
peak_exceed <= best_peak_exceed + 1e-3
total_exceed <= best_total_exceed + 1e-2

但这里的 exceed_* 是辅助变量，不是实际超限值本身。
它们只有：

python
exceed >= abs(x) - 0.5

而没有强制等于 max(abs(x)-0.5, 0)。

这意味着第二阶段可以这样干：

  • 把某些点真实弦值变差
  • 再把 exceed 在别的点上重新分配
  • 只要总和和峰值不超第一阶段上限，就还是“合法”

所以会出现一种情况：

  • 第一阶段确实把关键峰压下来了
  • 第二阶段为了顺坡，把别的位置又放坏了
  • 结果你看总体反而更差

这个问题是我那版两阶段写法里最像 bug 的地方。

2. 我把“0.5 一级目标”追得太狠了，但你的物理模型本身还没验证

这会导致一个典型现象：

原版模型虽然不够激进，但它至少是在一个“错误没被充分放大”的状态下工作。
我那版把 10m/20m 的 0.5 作为一级目标狠狠干，等于把模型误差也一起放大了。

如果下面任何一条不对，结果都会更差：

  • STEP_20 = 17 不是真实 20m 弦离散步长
  • 10m/20m 弦的公式跟现场检测系统不一致
  • 轨向/高低不该只由左股承担
  • 右股在现场其实也参与消峰，但模型里没参与

原版因为只是“尽量优化”，模型错一点还能勉强看。
我那版因为“死盯 0.5”，模型只要有一点偏，结果就会明显变差。

3. 你的原始模型本来就只有“左股在修弦值”，这件事在强目标下会更吃亏

你原始代码里这几句没变：

python
new_l10h = raw_10h - (A10 @ h_l)
new_l20h = raw_20h - (A20 @ h_l)

new_l10v = raw_10v + (A10 @ v_l)
new_l20v = raw_20v + (A20 @ v_l)

也就是说：

  • 轨向 10m/20m 全靠 h_l
  • 高低 10m/20m 全靠 v_l

右股只参与：

python
new_gauge = gauge_org - h_l + h_r
new_level = level_org - v_l + v_r

这等于告诉优化器：

  • 左股负责消波
  • 右股只负责配平轨距/水平

但现场人工往往不是这么干的，很多时候是左右股联动分摊。
所以我那版越强调“把弦值狠狠干到 0.5”，越会暴露这个自由度不足的问题。

4. 0.5 目标对“0.5mm 输出取整”非常敏感

这个点很容易被忽略。

你保存时是：

python
np.round(value * 2) / 2

也就是最终都按 0.5 mm 取整。

原版目标比较松，取整后恶化不一定太明显。
我那版很多点可能是靠 0.1、0.2、0.3 mm 的细微分配，才勉强把弦值压在 0.5 附近。
一旦最后统一取 0.5，可能直接把不少点又弹回去。

更关键的是：
你打印的报告是对“未取整解”做的，不是对导出后的取整方案做的。

所以有可能出现：

  • 求解器内部解看起来很好
  • 导出 CSV 后重新验算，明显变差

这个问题两版都有，但我那版会更敏感。

5. 第一阶段和第二阶段的目标之间有“风格冲突”

原版其实是单一风格：

  • 顺坡优先
  • 少动道优先
  • 指标尽量压

我那版变成：

  • 第一阶段：主峰优先，狠狠干超 0.5
  • 第二阶段：强顺坡、强头尾过渡、少动道

这两种风格不一定兼容。
如果第一阶段为了打掉某些峰，已经形成了比较“现场式”的局部动作，第二阶段可能又把它抹平。

所以它不是简单的“两阶段更高级”，而是可能变成：

第一阶段像老师傅，第二阶段像设计院，最后互相打架。

我现在最怀疑的结论

如果只给一个判断，我会说：

我那版变差，最可能是这三件事叠加：

  1. 第二阶段保护第一阶段成果的方式不严密
  2. 0.5 一级目标把原模型误差放大了
  3. 最终 0.5mm 取整把本来就脆弱的优化结果又破坏了

最该立刻做的 3 个定位实验

实验 1：分别看“第一阶段结束”和“第二阶段结束”的指标

如果：

  • 第一阶段好
  • 第二阶段差

那就基本坐实：第二阶段把第一阶段冲掉了。

实验 2：分别看“取整前”和“取整后”的指标

如果：

  • 取整前不错
  • 取整后明显恶化

那就说明：0.5mm 离散化是主因之一。

实验 3：只保留 10m 一级目标，先别把 20m 也狠狠干

如果：

  • 只追 10m 时效果明显变好
  • 加上 20m 后变差

那就说明：20m 模型或 20m 权重/步长有问题，尤其是 STEP_20 = 17 很可疑。

我会优先怎么修

不是再盲调权重，而是先做这三件实改：

第一，把第二阶段的“锁”改成锁真实结果，不锁辅助变量

也就是第二阶段直接约束：

python
cp.abs(new_l10h[valid_10]) <= stage1_best_10h + tol
cp.abs(new_l20h[valid_20]) <= stage1_best_20h + tol
...

或者锁住每类指标的最大值、均方值，而不是锁 exceed 变量。

第二，导出前加入“取整后复验”

把 round_to_0.5 后的 h_l/h_r/v_l/v_r 重新代回去算一次，别只看连续解。

第三，先做“1.5阶段”

也就是：

  • 第一阶段压 10m 到 0.5
  • 第二阶段只轻微顺坡
  • 暂时不要把 20m 也作为一级强目标

这样更稳。

直说

有头绪，而且我觉得问题不主要在“权重还要再调”，而在于：

我上次那版把两阶段写成了“形式上两阶段”，但没有把第一阶段成果真正严密地保护住。

这属于建模方式的问题，不是简单参数问题。

下一步最值钱的事不是继续猜，而是我直接帮你改成一个**“阶段1结果强锁定 + 取整后复验”的修正版**。

سؤال

改成一个**“阶段1结果强锁定 + 取整后复验”的修正版**。

الفكرة

我在修两处真正关键的地方：一是第二阶段直接锁“真实弦值结果”，不再锁辅助变量；二是导出前把 0.5 mm 取整后的方案重新代回模型复验，避免“连续解很好、落地方案变差”。

Thought for 46s

الإجابة

可以，下面这版就是你要的修正版，只做两件关键修正：

  1. 阶段 2 直接强锁定阶段 1 的真实结果
    不再锁 exceed 这种辅助变量，而是直接锁：

    • |10m轨向| 每个点
    • |20m轨向| 每个点
    • |10m高低| 每个点
    • |20m高低| 每个点

    这样第二阶段就不能把第一阶段已经压下来的点再放坏。

  2. 导出前做 0.5 mm 取整后复验
    连续解和最终写入 CSV 的离散解，分别都重新验算一遍。
    这样你看到的报告就是“真实落地方案”的效果，不会再出现“求解器里很好，导出后变差但你没看到”的情况。

直接替换版代码

python
import os
import sys

import cvxpy as cp
import numpy as np
import pandas as pd


# =========================
# 参数区
# =========================
TARGET_LIMIT = 0.5      # 一级目标：10m / 20m 尽量压到 0.5 mm 以内
BACKUP_LIMIT = 1.2      # 兜底硬限
GAUGE_LIMIT = 1.0       # 轨距偏差硬限
LEVEL_LIMIT = 1.2       # 水平偏差硬限
TRIANGLE_LIMIT = 1.5    # 三角坑硬限

STEP_10 = 8
STEP_20 = 17            # 这个必须与你的点距/检测公式核对

MAX_LIFT = 15.0
MAX_LINE = 5.0

MAX_DIFF_H = 0.6        # 相邻点拨道变化限值
MAX_DIFF_V = 0.6        # 相邻点起道变化限值

ROUND_PRECISION = 0.5   # 输出按 0.5 mm 取整

# 阶段2强锁定容差：非常小，但给一点数值余量
STRONG_LOCK_TOL = 0.02  # mm


def read_csv_auto(path):
    try:
        return pd.read_csv(path, encoding="gbk")
    except Exception:
        return pd.read_csv(path, encoding="utf-8")


def get_numeric_column(df, idx):
    """
    数值列清洗：
    1. 非数值转 NaN
    2. 插值补缺
    3. 首尾再补齐
    """
    s = pd.to_numeric(df.iloc[:, idx], errors="coerce").astype(float)
    if s.notna().sum() == 0:
        raise ValueError(f"第 {idx} 列没有有效数值。")
    s = s.interpolate(limit_direction="both")
    s = s.ffill().bfill().fillna(0.0)
    return s.to_numpy(dtype=float)


def build_matrix(size, step):
    """
    弦测法离散算子：
    y_i - 0.5 * (y_{i-step} + y_{i+step})
    """
    main_diag = np.ones(size)
    side_diag = -0.5 * np.ones(size - step)
    return (
        np.diag(main_diag)
        + np.diag(side_diag, k=step)
        + np.diag(side_diag, k=-step)
    )


def round_to_step(x, step=0.5):
    return np.round(x / step) * step


def solve_problem(prob, stage_name):
    """
    优先 OSQP，失败时回退 SCS
    """
    tried = []

    try:
        prob.solve(
            solver=cp.OSQP,
            verbose=False,
            warm_start=True,
            max_iter=100000,
            eps_abs=1e-4,
            eps_rel=1e-4,
        )
        tried.append(("OSQP", prob.status))
        if prob.status in ("optimal", "optimal_inaccurate"):
            print(f"{stage_name} 求解成功: OSQP / {prob.status}")
            return
    except Exception as e:
        tried.append(("OSQP", str(e)))

    try:
        prob.solve(
            solver=cp.SCS,
            verbose=False,
            warm_start=True,
            max_iters=40000,
            eps=1e-4,
        )
        tried.append(("SCS", prob.status))
        if prob.status in ("optimal", "optimal_inaccurate"):
            print(f"{stage_name} 求解成功: SCS / {prob.status}")
            return
    except Exception as e:
        tried.append(("SCS", str(e)))

    raise RuntimeError(f"{stage_name} 求解失败，尝试记录: {tried}")


def evaluate_scheme(
    raw_10h,
    raw_20h,
    gauge_org,
    raw_10v,
    raw_20v,
    level_org,
    h_l,
    h_r,
    v_l,
    v_r,
    A10,
    A20,
    valid_10,
    valid_20,
):
    """
    对一个具体方案（连续解或取整解）做复验
    """
    new_gauge = gauge_org - h_l + h_r
    new_level = level_org - v_l + v_r

    new_l10h = raw_10h - (A10 @ h_l)
    new_l20h = raw_20h - (A20 @ h_l)

    new_l10v = raw_10v + (A10 @ v_l)
    new_l20v = raw_20v + (A20 @ v_l)

    abs_10h = np.abs(new_l10h[valid_10])
    abs_20h = np.abs(new_l20h[valid_20])
    abs_10v = np.abs(new_l10v[valid_10])
    abs_20v = np.abs(new_l20v[valid_20])

    triangle = np.abs(new_level[5:] - new_level[:-5])

    def ratio_leq(arr, th):
        return 100.0 * np.mean(arr <= th)

    out = {
        "new_gauge": new_gauge,
        "new_level": new_level,
        "new_l10h": new_l10h,
        "new_l20h": new_l20h,
        "new_l10v": new_l10v,
        "new_l20v": new_l20v,
        "abs_10h": abs_10h,
        "abs_20h": abs_20h,
        "abs_10v": abs_10v,
        "abs_20v": abs_20v,
        "triangle": triangle,
        "max_gauge": float(np.max(np.abs(new_gauge))),
        "max_level": float(np.max(np.abs(new_level))),
        "max_10h": float(np.max(abs_10h)),
        "max_20h": float(np.max(abs_20h)),
        "max_10v": float(np.max(abs_10v)),
        "max_20v": float(np.max(abs_20v)),
        "max_triangle": float(np.max(triangle)),
        "pct_10h_le_05": ratio_leq(abs_10h, 0.5),
        "pct_20h_le_05": ratio_leq(abs_20h, 0.5),
        "pct_10v_le_05": ratio_leq(abs_10v, 0.5),
        "pct_20v_le_05": ratio_leq(abs_20v, 0.5),
        "tqi_10h_before": float(np.std(raw_10h[valid_10])),
        "tqi_20h_before": float(np.std(raw_20h[valid_20])),
        "tqi_10v_before": float(np.std(raw_10v[valid_10])),
        "tqi_20v_before": float(np.std(raw_20v[valid_20])),
        "tqi_10h_after": float(np.std(new_l10h[valid_10])),
        "tqi_20h_after": float(np.std(new_l20h[valid_20])),
        "tqi_10v_after": float(np.std(new_l10v[valid_10])),
        "tqi_20v_after": float(np.std(new_l20v[valid_20])),
    }
    return out


def print_report(title, rep):
    print("\n" + "=" * 62)
    print(title)
    print(f"最大轨距偏差: {rep['max_gauge']:.3f} mm")
    print(f"最大水平偏差: {rep['max_level']:.3f} mm")
    print(f"最大10m轨向:   {rep['max_10h']:.3f} mm")
    print(f"最大20m轨向:   {rep['max_20h']:.3f} mm")
    print(f"最大10m高低:   {rep['max_10v']:.3f} mm")
    print(f"最大20m高低:   {rep['max_20v']:.3f} mm")
    print(f"最大三角坑:    {rep['max_triangle']:.3f} mm")
    print("-" * 62)
    print(f"10m轨向 <= 0.5 占比: {rep['pct_10h_le_05']:6.2f}%")
    print(f"20m轨向 <= 0.5 占比: {rep['pct_20h_le_05']:6.2f}%")
    print(f"10m高低 <= 0.5 占比: {rep['pct_10v_le_05']:6.2f}%")
    print(f"20m高低 <= 0.5 占比: {rep['pct_20v_le_05']:6.2f}%")
    print("-" * 62)
    print(f"10m轨向 TQI: {rep['tqi_10h_before']:.3f} -> {rep['tqi_10h_after']:.3f}")
    print(f"20m轨向 TQI: {rep['tqi_20h_before']:.3f} -> {rep['tqi_20h_after']:.3f}")
    print(f"10m高低 TQI: {rep['tqi_10v_before']:.3f} -> {rep['tqi_10v_after']:.3f}")
    print(f"20m高低 TQI: {rep['tqi_20v_before']:.3f} -> {rep['tqi_20v_after']:.3f}")
    print("=" * 62)


def run_two_stage_manual_style_tuning_strong_lock(input_path, output_path):
    # =========================
    # 1. 数据读取
    # =========================
    df = read_csv_auto(input_path)

    raw_10h = get_numeric_column(df, 3)    # 原始10m轨向
    raw_20h = get_numeric_column(df, 5)    # 原始20m轨向
    gauge_org = get_numeric_column(df, 7)  # 原始轨距偏差

    raw_10v = get_numeric_column(df, 8)    # 原始10m高低
    raw_20v = get_numeric_column(df, 10)   # 原始20m高低
    level_org = get_numeric_column(df, 12) # 原始水平偏差

    n = len(df)
    if n <= 2 * STEP_20:
        raise ValueError(f"数据点数不足，当前 {n} 点，无法形成 20m 有效区。")

    valid_10 = slice(STEP_10, -STEP_10)
    valid_20 = slice(STEP_20, -STEP_20)

    m10 = n - 2 * STEP_10
    m20 = n - 2 * STEP_20

    A10 = build_matrix(n, STEP_10)
    A20 = build_matrix(n, STEP_20)

    # =========================
    # 2. 决策变量
    # =========================
    v_l = cp.Variable(n)
    v_r = cp.Variable(n)
    h_l = cp.Variable(n)
    h_r = cp.Variable(n)

    # =========================
    # 3. 物理模型
    # =========================
    new_gauge = gauge_org - h_l + h_r
    new_level = level_org - v_l + v_r

    new_l10h = raw_10h - (A10 @ h_l)
    new_l20h = raw_20h - (A20 @ h_l)

    new_l10v = raw_10v + (A10 @ v_l)
    new_l20v = raw_20v + (A20 @ v_l)

    # =========================
    # 4. 基础约束
    # =========================
    base_constraints = [
        # 量程
        v_l >= 0,
        v_l <= MAX_LIFT,
        v_r >= 0,
        v_r <= MAX_LIFT,
        h_l >= -MAX_LINE,
        h_l <= MAX_LINE,
        h_r >= -MAX_LINE,
        h_r <= MAX_LINE,

        # 安全硬约束
        new_gauge >= -GAUGE_LIMIT,
        new_gauge <= GAUGE_LIMIT,
        new_level >= -LEVEL_LIMIT,
        new_level <= LEVEL_LIMIT,

        # 三角坑
        cp.abs(new_level[5:] - new_level[:-5]) <= TRIANGLE_LIMIT,

        # 相邻点变化
        cp.abs(cp.diff(h_l)) <= MAX_DIFF_H,
        cp.abs(cp.diff(h_r)) <= MAX_DIFF_H,
        cp.abs(cp.diff(v_l)) <= MAX_DIFF_V,
        cp.abs(cp.diff(v_r)) <= MAX_DIFF_V,
    ]

    # =========================
    # 5. 第一阶段：先压 0.5 目标
    # =========================
    exceed_10h = cp.Variable(m10, nonneg=True)
    exceed_10v = cp.Variable(m10, nonneg=True)
    exceed_20h = cp.Variable(m20, nonneg=True)
    exceed_20v = cp.Variable(m20, nonneg=True)

    backup_10h = cp.Variable(m10, nonneg=True)
    backup_10v = cp.Variable(m10, nonneg=True)
    backup_20h = cp.Variable(m20, nonneg=True)
    backup_20v = cp.Variable(m20, nonneg=True)

    peak_exceed = cp.Variable(nonneg=True)

    stage1_constraints = base_constraints + [
        # exceed >= |x| - 0.5
        exceed_10h >= cp.abs(new_l10h[valid_10]) - TARGET_LIMIT,
        exceed_10v >= cp.abs(new_l10v[valid_10]) - TARGET_LIMIT,
        exceed_20h >= cp.abs(new_l20h[valid_20]) - TARGET_LIMIT,
        exceed_20v >= cp.abs(new_l20v[valid_20]) - TARGET_LIMIT,

        # 主峰优先
        exceed_10h <= peak_exceed,
        exceed_10v <= peak_exceed,
        exceed_20h <= peak_exceed,
        exceed_20v <= peak_exceed,

        # 兜底硬限
        cp.abs(new_l10h[valid_10]) <= BACKUP_LIMIT + backup_10h,
        cp.abs(new_l10v[valid_10]) <= BACKUP_LIMIT + backup_10v,
        cp.abs(new_l20h[valid_20]) <= BACKUP_LIMIT + backup_20h,
        cp.abs(new_l20v[valid_20]) <= BACKUP_LIMIT + backup_20v,
    ]

    stage1_obj = cp.Minimize(
        1_000_000.0 * peak_exceed
        + 20_000.0 * (cp.sum(exceed_10h) + cp.sum(exceed_10v))
        + 12_000.0 * (cp.sum(exceed_20h) + cp.sum(exceed_20v))
        + 300_000.0 * (
            cp.sum(backup_10h) + cp.sum(backup_10v)
            + cp.sum(backup_20h) + cp.sum(backup_20v)
        )
        # 只保留很轻的正则，防止第一阶段解发散
        + 10.0 * cp.sum_squares(h_l)
        + 10.0 * cp.sum_squares(h_r)
        + 4.0 * (cp.sum(v_l) + cp.sum(v_r))
        + 20.0 * cp.sum_squares(cp.diff(h_l, 2))
        + 20.0 * cp.sum_squares(cp.diff(v_l, 2))
    )

    print(f"正在执行第一阶段：先压 10m/20m 到 0.5 附近（点数 {n}）...")
    prob1 = cp.Problem(stage1_obj, stage1_constraints)
    solve_problem(prob1, "第一阶段")

    if h_l.value is None:
        raise RuntimeError("第一阶段未找到可行解。")

    # 第一阶段真实结果，用于第二阶段强锁定
    stage1_rep = evaluate_scheme(
        raw_10h, raw_20h, gauge_org,
        raw_10v, raw_20v, level_org,
        h_l.value, h_r.value, v_l.value, v_r.value,
        A10, A20, valid_10, valid_20
    )

    print_report("第一阶段复验报告（连续解）", stage1_rep)

    # 把第一阶段每个点的真实绝对值锁下来
    lock_10h = stage1_rep["abs_10h"]
    lock_20h = stage1_rep["abs_20h"]
    lock_10v = stage1_rep["abs_10v"]
    lock_20v = stage1_rep["abs_20v"]

    # =========================
    # 6. 第二阶段：强锁定阶段1结果，再顺坡/少动道
    # =========================
    stage2_constraints = base_constraints + [
        # 关键修正：直接锁真实结果，不锁辅助变量
        cp.abs(new_l10h[valid_10]) <= lock_10h + STRONG_LOCK_TOL,
        cp.abs(new_l20h[valid_20]) <= lock_20h + STRONG_LOCK_TOL,
        cp.abs(new_l10v[valid_10]) <= lock_10v + STRONG_LOCK_TOL,
        cp.abs(new_l20v[valid_20]) <= lock_20v + STRONG_LOCK_TOL,
    ]

    stage2_obj = cp.Minimize(
        # A. 顺坡主目标
        3000.0 * cp.sum_squares(cp.diff(h_l, 2))
        + 3000.0 * cp.sum_squares(cp.diff(v_l, 2))
        + 900.0 * cp.sum_squares(cp.diff(h_r, 2))
        + 900.0 * cp.sum_squares(cp.diff(v_r, 2))

        # B. 轨距 / 水平优化
        + 700.0 * cp.sum_squares(new_gauge)
        + 700.0 * cp.sum_squares(new_level)

        # C. 少动道
        + 18.0 * cp.sum_squares(h_l)
        + 18.0 * cp.sum_squares(h_r)
        + 5.0 * (cp.sum(v_l) + cp.sum(v_r))

        # D. 头尾过渡
        + 180.0 * (
            cp.sum_squares(h_l[:STEP_20]) + cp.sum_squares(h_l[-STEP_20:])
            + cp.sum_squares(h_r[:STEP_20]) + cp.sum_squares(h_r[-STEP_20:])
            + cp.sum_squares(v_l[:STEP_20]) + cp.sum_squares(v_l[-STEP_20:])
            + cp.sum_squares(v_r[:STEP_20]) + cp.sum_squares(v_r[-STEP_20:])
        )
    )

    print("正在执行第二阶段：强锁定阶段1结果，再做顺坡/少动道...")
    prob2 = cp.Problem(stage2_obj, stage2_constraints)
    solve_problem(prob2, "第二阶段")

    if h_l.value is None:
        raise RuntimeError("第二阶段未找到可行解。")

    # =========================
    # 7. 连续解复验
    # =========================
    cont_rep = evaluate_scheme(
        raw_10h, raw_20h, gauge_org,
        raw_10v, raw_20v, level_org,
        h_l.value, h_r.value, v_l.value, v_r.value,
        A10, A20, valid_10, valid_20
    )

    print_report("第二阶段复验报告（连续解）", cont_rep)

    # =========================
    # 8. 取整后复验
    # =========================
    h_l_round = round_to_step(h_l.value, ROUND_PRECISION)
    h_r_round = round_to_step(h_r.value, ROUND_PRECISION)
    v_l_round = round_to_step(v_l.value, ROUND_PRECISION)
    v_r_round = round_to_step(v_r.value, ROUND_PRECISION)

    round_rep = evaluate_scheme(
        raw_10h, raw_20h, gauge_org,
        raw_10v, raw_20v, level_org,
        h_l_round, h_r_round, v_l_round, v_r_round,
        A10, A20, valid_10, valid_20
    )

    print_report("最终导出方案复验报告（0.5mm 取整后）", round_rep)

    # 取整恶化提示
    print("\n取整影响检查：")
    print(f"10m轨向最大值: 连续解 {cont_rep['max_10h']:.3f} -> 取整后 {round_rep['max_10h']:.3f}")
    print(f"20m轨向最大值: 连续解 {cont_rep['max_20h']:.3f} -> 取整后 {round_rep['max_20h']:.3f}")
    print(f"10m高低最大值: 连续解 {cont_rep['max_10v']:.3f} -> 取整后 {round_rep['max_10v']:.3f}")
    print(f"20m高低最大值: 连续解 {cont_rep['max_20v']:.3f} -> 取整后 {round_rep['max_20v']:.3f}")

    if (
        round_rep["max_gauge"] > GAUGE_LIMIT + 1e-6
        or round_rep["max_level"] > LEVEL_LIMIT + 1e-6
        or round_rep["max_triangle"] > TRIANGLE_LIMIT + 1e-6
    ):
        print("[警告] 取整后，安全类指标出现越限，请重点检查导出方案。")

    if (
        round_rep["max_10h"] > cont_rep["max_10h"] + 0.2
        or round_rep["max_20h"] > cont_rep["max_20h"] + 0.2
        or round_rep["max_10v"] > cont_rep["max_10v"] + 0.2
        or round_rep["max_20v"] > cont_rep["max_20v"] + 0.2
    ):
        print("[提示] 取整对弦值破坏较明显，说明方案对 0.5mm 离散化比较敏感。")

    # =========================
    # 9. 写回 CSV（写入的是最终取整方案）
    # =========================
    df.iloc[:, 13] = h_l_round
    df.iloc[:, 14] = h_r_round
    df.iloc[:, 15] = v_l_round
    df.iloc[:, 16] = v_r_round

    df.to_csv(output_path, index=False, encoding="gbk")


if __name__ == "__main__":
    input_file = ""

    if len(sys.argv) > 1:
        input_file = sys.argv[1]
    else:
        print("--- 铁路轨道两阶段精调工具（强锁定 + 取整复验版）---")
        print("提示：你可以直接将 CSV 文件拖入此窗口，然后按回车")
        input_file = input("请输入或拖入要处理的 CSV 文件路径: ").strip()
        input_file = input_file.replace('"', "").replace("'", "")

    if input_file:
        if os.path.exists(input_file) and input_file.lower().endswith(".csv"):
            base, ext = os.path.splitext(input_file)
            output_file = f"{base}_两阶段强锁定优化方案{ext}"

            try:
                run_two_stage_manual_style_tuning_strong_lock(input_file, output_file)
                print(f"\n[成功] 方案已生成：{os.path.basename(output_file)}")
            except Exception as e:
                print(f"\n[错误] 运行失败，详细信息: {e}")
        else:
            print(f"\n[错误] 找不到文件或文件格式不是 CSV: {input_file}")
    else:
        print("\n[错误] 未输入任何有效路径。")

    print("\n" + "-" * 45)
    input("按回车键退出程序...")

这版和上一版相比，真正修了什么

1. 阶段 2 的锁定方式变了

上一版的问题是锁了辅助变量，第二阶段仍有机会“换地方变差”。

现在改成直接锁真实结果：

python
cp.abs(new_l10h[valid_10]) <= lock_10h + STRONG_LOCK_TOL
cp.abs(new_l20h[valid_20]) <= lock_20h + STRONG_LOCK_TOL
cp.abs(new_l10v[valid_10]) <= lock_10v + STRONG_LOCK_TOL
cp.abs(new_l20v[valid_20]) <= lock_20v + STRONG_LOCK_TOL

这意味着：

  • 第一阶段某个点是 0.42
  • 第二阶段最多只能变成 0.44 左右
  • 不可能再被顺坡目标“冲坏”

这才叫真正的强锁定。

2. 现在会打印三份报告

你会同时看到：

  • 第一阶段复验报告（连续解）
  • 第二阶段复验报告（连续解）
  • 最终导出方案复验报告（0.5mm 取整后）

这样你可以立刻定位：

  • 是第一阶段本来就没压好
  • 还是第二阶段冲坏了
  • 还是取整以后坏掉了

你跑完以后，先看哪三行最有价值

先看这三组对比：

第一组

第一阶段复验报告（连续解）

如果这里已经很差，说明：

  • 不是第二阶段问题
  • 而是一级目标本身就没打下来
  • 下一步该怀疑物理模型或 20m 步长

第二组

第二阶段复验报告（连续解）

如果第一阶段好、第二阶段仍几乎一样，那说明：

  • 强锁定生效了
  • 第二阶段没有再冲坏

第三组

最终导出方案复验报告（0.5mm 取整后）

如果这里明显比连续解差，说明：

  • 不是优化器的问题
  • 是0.5 mm 离散化本身在破坏方案

我先提前告诉你一个预期

这版如果仍然比你最初原版差，那最值得怀疑的就不再是“两阶段写法”了，而是下面两个结构性问题：

  1. STEP_20 = 17 不是你的真实 20m 弦离散步长
  2. 你的 10m/20m 指标并不是“只由左股起拨道单独决定”

也就是说，到那一步，问题就基本不在“优化技巧”，而在物理映射模型本身了。

你把这版跑一下，然后把这三份报告里这几项贴给我就够了：

  • 第一阶段连续解的 最大10m轨向 / 最大20m轨向 / 最大10m高低 / 最大20m高低
  • 第二阶段连续解的这四项
  • 取整后的这四项

我就能很快判断，问题到底卡在“阶段1”“阶段2”还是“取整”。

مشاركة هذا السؤال والجواب