如何将求解仪器时钟漂移校正方程组的SymPy代码转为SciPy版本

我需要求解如下形式的方程组：

其中已知量为c、a_ir以及向量p_r，未知量为向量x0、v以及所有ti参数。另外T和R为已知但会动态变化，无法硬编码。原实现采用SymPy的nsolve性能不足，需要转换为SciPy兼容的形式完成求解。

SciPy的非线性求解器依赖残差函数作为输入：输入为所有未知量拼接的一维数组，输出为所有方程的残差（等式左侧减右侧，求解目标为残差全为0）。可以直接使用SymPy的lambdify接口把已有的符号表达式自动转换为数值函数，不需要手动改写公式避免出错。

import numpy as np
import sympy as sp
from dataclasses import dataclass
from scipy.optimize import root

# 按实际需求替换为你的c值
c = 3e8

@dataclass
class Vertexer:
    receivers: list

    def __post_init__(self):
        self.R = len(self.receivers)
        # 固定符号变量提前定义，避免每次求解重复生成
        self.x, self.y, self.z, self.vx, self.vy, self.vz = sp.symbols('x y z vx vy vz', real=True)
        self.base_vars = (self.x, self.y, self.z, self.vx, self.vy, self.vz)

    def find(self, data, init_guess=None):
        R = self.R
        T = len(data[0])
        # 动态生成ti符号变量
        t_syms = sp.symbols(f't0:{T}', real=True)
        all_syms = self.base_vars + t_syms
        # 生成所有方程的残差（左减右）
        residuals = []
        for r in range(R):
            prx, pry, prz = self.receivers[r]
            dx = self.x - prx
            dy = self.y - pry
            dz = self.z - prz
            v_sq = self.vx**2 + self.vy**2 + self.vz**2
            v_dot_dx = self.vx * dx + self.vy * dy + self.vz * dz
            for i in range(T):
                a_ir = data[r][i]
                ti = t_syms[i]
                t0 = t_syms[0]
                delta_t = ti - t0
                left = c**2 * (a_ir - ti)**2
                right = v_sq * delta_t**2 + 2 * delta_t * v_dot_dx + (dx**2 + dy**2 + dz**2)
                residuals.append(left - right)
        # 转换为数值计算函数
        f_res = sp.lambdify(all_syms, residuals, 'numpy')
        # 封装为SciPy要求的输入格式：输入一维数组，输出一维残差数组
        def objective(x):
            return f_res(*x)
        # 未传入初始猜测时生成默认初值，可根据业务逻辑调整
        if init_guess is None:
            # x0初始值设为所有接收站的平均位置
            mean_rx = np.mean([p[0] for p in self.receivers])
            mean_ry = np.mean([p[1] for p in self.receivers])
            mean_rz = np.mean([p[2] for p in self.receivers])
            # v初始值可按实际场景调整，示例为光速的一半
            init_v = c/2 * np.ones(3)
            # ti初始值设为所有a_ir对应位置的平均值
            init_t = np.mean(data, axis=0)
            init_guess = np.concatenate([[mean_rx, mean_ry, mean_rz], init_v, init_t])
        # 调用求解器
        sol = root(objective, init_guess, method='lm')
        # 拆分求解结果
        x0 = sol.x[:3]
        v = sol.x[3:6]
        ti = sol.x[6:]
        return {
            'success': sol.success,
            'x0': x0,
            'v': v,
            'ti': ti,
            'residual_norm': np.linalg.norm(sol.fun)
        }

• 初始化时传入接收站坐标列表receivers，格式为[[p1x,p1y,p1z], [p2x,p2y,p2z], ...]
• 调用find方法时传入data，格式为R行T列的二维数组，对应所有a_ir值
• 可选传入init_guess参数自定义初始猜测值，非线性求解对初值敏感，合理的初值可以大幅提升求解成功率和速度
• 求解结果中可通过success字段判断是否收敛，其余字段分别对应各未知量的求解结果和残差二范数

内容的提问来源于stack exchange，提问作者10GeV