使用GeoSolver与Sympy求解3D约束时如何仅获取实数解？

问题背景 #

我用GeoSolver包求解3D约束，处理平行约束时得到了复数解，但实际存在实数解，想知道怎么只获取实数解。

当前求解代码 #

result = sp.solve(AllEquestions, AllVariables)

约束方程与变量 #

• 约束方程AllEquestions：

['x1- 0', 'y1- 0', 'z1- 0', ((-x1 + x2)*(-x3 + x4) + (-y1 + y2)*(-y3 + y4) + (-z1 + z2)*(-z3 + z4))**2/(((x1 - x2)**2 + (y1 - y2)**2 + (z1 - z2)**2)**1.0*((x3 - x4)**2 + (y3 - y4)**2 + (z3 - z4)**2)**1.0) - 1, 'x2- 100', 'y2- 25', 'z2- 0', ((-x1 + x2)*(-x3 + x4) + (-y1 + y2)*(-y3 + y4) + (-z1 + z2)*(-z3 + z4))**2/(((x1 - x2)**2 + (y1 - y2)**2 + (z1 - z2)**2)**1.0*((x3 - x4)**2 + (y3 - y4)**2 + (z3 - z4)**2)**1.0) - 1, 'x3- 10', 'y3- 10', 'z3- 0', ((-x1 + x2)*(-x3 + x4) + (-y1 + y2)*(-y3 + y4) + (-z1 + z2)*(-z3 + z4))**2/(((x1 - x2)**2 + (y1 - y2)**2 + (z1 - z2)**2)**1.0*((x3 - x4)**2 + (y3 - y4)**2 + (z3 - z4)**2)**1.0) - 1, 'x4- 35', ((-x1 + x2)*(-x3 + x4) + (-y1 + y2)*(-y3 + y4) + (-z1 + z2)*(-z3 + z4))**2/(((x1 - x2)**2 + (y1 - y2)**2 + (z1 - z2)**2)**1.0*((x3 - x4)**2 + (y3 - y4)**2 + (z3 - z4)**2)**1.0) - 1]

• 变量AllVariables：

[x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4]

当前返回结果 #

[(0.0, 0.0, 0.0, 100.000000000000, 25.0000000000000, 0.0, 10.0000000000000, 10.0000000000000, 0.0, 35.0000000000000, -1.03077640640442*I*z4 + 16.25, z4), (0.0, 0.0, 0.0, 100.000000000000, 25.0000000000000, 0.0, 10.0000000000000, 10.0000000000000, 0.0, 35.0000000000000, 1.03077640640442*I*z4 + 16.25, z4)]

期望结果 #

[(0.0, 0.0, 0.0, 100.000000000000, 25.0000000000000, 0.0, 10.0000000000000, 10.0000000000000, 0.0, 35.0000000000000, 16.2500000000000, 0.0)]

尝试过的方法及错误 #

原本定义变量的方式：

self.x = sp.Symbol('x' + str(self.local_var))

为变量添加real=True标志后：

self.x = sp.Symbol('x' + str(self.local_var), real=True)

触发错误：

File "C:\Users\Achaibou Karim\AppData\Roaming\FreeCAD\Macro\fpo\tube\GeoSolver\solver.py", line 25, in solve
result = sp.solve(AllEquestions, AllVariables)
File "c:\Program Files\FreeCAD 0.21\bin\lib\site-packages\sympy\solvers\solvers.py", line 1172, in solve
linear, solution = _solve_system(f, symbols, **flags)
File "c:\Program Files\FreeCAD 0.21\bin\lib\site-packages\sympy\solvers\solvers.py", line 1896, in _solve_system
raise NotImplementedError('no valid subset found')
NotImplementedError: no valid subset found

解决方案 #

方法1：后处理过滤实数解 #

不需要修改变量定义，在得到解之后手动过滤掉包含复数的结果，结合你的期望结果，可直接指定z4=0来获取目标实数解：

real_results = []
for sol in result:
    z4_val = sol[11]
    y4_val = sol[10]
    # 当z4=0时，y4为实数16.25，符合期望
    if z4_val == 0 and y4_val.is_real:
        real_results.append(sol)

# 通用过滤逻辑：检查解中所有变量是否为实数
# real_results = [sol for sol in result if all(isinstance(v, (int, float)) or v.is_real for v in sol)]

方法2：优化平行约束方程 #

当前用的点积平方分式方程等价于向量平行，但浮点数指数**1.0会增加Sympy的处理难度。可以替换为向量叉乘为零向量的形式，这是更直接的多项式约束，Sympy更容易处理实数解：
原平行约束替换为以下三个方程（对应叉乘的三个分量）：

# 向量1：(x2-x1, y2-y1, z2-z1)
# 向量2：(x4-x3, y4-y3, z4-z3)
(y2 - y1)*(z4 - z3) - (z2 - z1)*(y4 - y3),
(z2 - z1)*(x4 - x3) - (x2 - x1)*(z4 - z3),
(x2 - x1)*(y4 - y3) - (y2 - y1)*(x4 - x3)

用这三个方程替换原来的分式约束，再配合real=True定义变量，Sympy应该能直接输出实数解，避免复数解生成。

方法3：调用solve时指定实数解参数 #

在调用solve函数时添加real=True参数，强制求解实数解：

result = sp.solve(AllEquestions, AllVariables, real=True)

注意：该方法对复杂分式方程的支持有限，优先推荐优化约束方程的方式。

内容的提问来源于stack exchange，提问作者Achaibou Karim