使用SymPy求解三角函数方程组返回EmptySet问题求助

为实现Pro-Link摩托车优化，我编写的Matlab代码通过vpasolve能快速求解三角函数方程组，但迁移到Python用SymPy库实现相同逻辑时，运行耗时约30秒且返回EmptySet，需要排查解决。

Matlab代码 #

clc
clear 

Lmono= 320;
Lbielletta=145;
IPS= 16.03;
Pivot= -20;

R1x=(-46.72)-((Pivot)*sind(IPS));
R1y=((180.26)-((Pivot)*cosd(IPS)));
R_1x=(-43.52)-((Pivot)*sind(IPS));
R_1y=((-151.37)-((Pivot)*cosd(IPS)));
R4=60; 
R3=203.727;
Ip=100.1;
eta=36.79;


syms phi o 
eqns = [(Lbielletta)^2 == ((R3*cosd(phi)-R4*sind(o)+R_1x)^2)+((-R3*sind(phi)-R4*cosd(o)-R_1y)^2)

(Lmono)^2 == ((((R3*cosd(phi))-(R4*sind(o))-((Ip*cosd(eta+o)))+(R1x))^2) + (((R3*sind(phi))+(R4*cosd(o))-((Ip*sind(eta+o)))+(R1y))^2))];

[phi ,o]=vpasolve(eqns,[phi o]);

Python代码（原问题版本） #

import math as m

Lb = 145.0
Lm = 320.0
IPS= 16.03;
Pivot= -20.0;
R1x=(-46.72)-((Pivot)*m.sin(IPS))
R1y=((180.26)-((Pivot)*m.cos(IPS)))
R_1x=(-43.52)-((Pivot)*m.sin(IPS))
R_1y=((-151.37)-((Pivot)*m.cos(IPS)))
R4=60.0
R3=203.727
Ip=100.1
eta=36.79

import sympy as sym
from sympy import sin, cos
    
sym.init_printing()
phi,o = sym.symbols('phi,o')
f = sym.Eq(((R3*cos(phi)-R4*sin(o)+R_1x)**2)+((-R3*sin(phi)-R4*cos(o)-R_1y)**2),Lb**2)
g = sym.Eq(((((R3*cos(phi))-(R4*sin(o))-((Ip*cos(eta+o)))+(R1x))**2) + (((R3*sin(phi))+(R4*cos(o))-((Ip*sin(eta+o)))+(R1y))**2)),Lm**2)
print(sym.nonlinsolve([f,g],(phi,o)))

Python运行结果 #

runfile('C:/Users/Administrator/.spyder-py3/temp.py', wdir='C:/Users/Administrator/.spyder-py3')
EmptySet

问题排查与解决方案 #

核心问题 #

1. 角度单位不匹配：Matlab的sind/cosd是角度制计算，而Python的math.sin/math.cos、SymPy的sin/cos默认是弧度制。原Python代码直接用角度值传入三角函数，导致R1x、R1y等参数计算错误，方程组实际无解，因此返回EmptySet。
2. 求解函数选择错误：SymPy的nonlinsolve主打符号解，对于这类非线性数值方程组，应该用和Matlabvpasolve对应的nsolve（数值求解函数），效率更高且适合此类场景。

修正后的Python代码 #

import sympy as sym
from sympy import sin, cos, nsolve, deg2rad, rad2deg

# 将所有角度参数转换为弧度（对齐Matlab sind/cosd的逻辑）
IPS = deg2rad(16.03)
Pivot = -20.0
eta = deg2rad(36.79)

# 用SymPy三角函数计算参数（保证精度，也可转弧度后用math计算数值）
R1x = (-46.72) - (Pivot * sin(IPS))
R1y = (180.26) - (Pivot * cos(IPS))
R_1x = (-43.52) - (Pivot * sin(IPS))
R_1y = (-151.37) - (Pivot * cos(IPS))

# 常量参数
Lb = 145.0
Lm = 320.0
R4 = 60.0
R3 = 203.727
Ip = 100.1

# 定义符号变量
phi, o = sym.symbols('phi,o')

# 构建方程组（所有三角函数均为弧度制）
eq1 = sym.Eq(
    (R3 * cos(phi) - R4 * sin(o) + R_1x)**2 + (-R3 * sin(phi) - R4 * cos(o) - R_1y)**2,
    Lb**2
)
eq2 = sym.Eq(
    (R3 * cos(phi) - R4 * sin(o) - Ip * cos(eta + o) + R1x)**2 +
    (R3 * sin(phi) + R4 * cos(o) - Ip * sin(eta + o) + R1y)**2,
    Lm**2
)

# 使用nsolve数值求解，提供初始猜测值（可根据场景调整）
solution = nsolve([eq1, eq2], (phi, o), (0, 0))

# 输出结果，可选择弧度或角度格式
print(f"phi（弧度）: {solution[0].evalf()}")
print(f"phi（角度）: {rad2deg(solution[0]).evalf()}")
print(f"o（弧度）: {solution[1].evalf()}")
print(f"o（角度）: {rad2deg(solution[1]).evalf()}")

说明 #

• 角度转换：用deg2rad将IPS、eta转为弧度，确保和SymPy三角函数的单位一致，Matlab的sind本质是sin(deg2rad(x))，这里完全对齐逻辑。
• 求解函数：nsolve需要初始猜测值，示例用(0,0)，如果求解失败可尝试其他合理初始值（比如Matlab得到的结果）。
• 效率：nsolve的求解速度和Matlabvpasolve接近，不会出现30秒的耗时问题。

内容的提问来源于stack exchange，提问作者Sergio Azzolina