如何在Python中用两个带约束的圆弧帽拟合数据点

1. 隐式方程拟合避免数据丢失：不用显式的y(x)表达式，直接用圆弧的隐式方程计算残差。对每个数据点(x,y)，残差定义为x² + (y - c)² - R²，最小化所有残差的平方和，这样不管y的取值（只要满足y<0）都能覆盖，不会因为开根号丢失数据。
2. 约束条件化简与融入：先把给定的约束条件化简，减少独立变量数量：
   • 代入cos(sigma) = -c_beta/R_beta到约束方程，最终可化简为：R_alpha² = c_alpha² + R_beta² - c_beta²
   • 由此解出c_alpha = sqrt(R_alpha² - R_beta² + c_beta²)（根据数据特征确定正负根，这里假设c_alpha为正）
   • 这样独立变量只剩R_beta、c_beta、R_alpha三个，直接将约束融入目标函数，无需额外的约束优化参数。

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

# ---------------------- 1. 替换为你的实际数据 ----------------------
# 示例：模拟灰色、绿色数据（用户需替换成自己的datagrey、datagreen）
# 灰色数据：圆心(0, 4)，半径5，y<0的部分
theta_grey = np.linspace(np.pi + np.arcsin(4/5), 2*np.pi - np.arcsin(4/5), 100)
x_grey = 5 * np.sin(theta_grey)
y_grey = 4 + 5 * np.cos(theta_grey)
datagrey = np.column_stack((x_grey, y_grey)) + np.random.normal(0, 0.1, (100, 2))

# 绿色数据：圆心(0, -3)，半径4，y<0的部分
theta_green = np.linspace(0, np.pi, 100)
x_green = 4 * np.sin(theta_green)
y_green = -3 + 4 * np.cos(theta_green)
datagreen = np.column_stack((x_green, y_green)) + np.random.normal(0, 0.1, (100, 2))

# ---------------------- 2. 定义拟合目标函数 ----------------------
def objective(params):
    R_beta, c_beta, R_alpha = params
    
    # 根据约束计算c_alpha
    c_alpha_sq = R_alpha**2 - R_beta**2 + c_beta**2
    if c_alpha_sq < 0:
        return np.inf  # 不满足约束，返回无穷大终止当前迭代
    c_alpha = np.sqrt(c_alpha_sq)
    
    # 计算灰色数据残差平方和
    x_g, y_g = datagrey[:, 0], datagrey[:, 1]
    res_grey = (x_g**2 + (y_g - c_alpha)**2 - R_alpha**2)**2
    sum_grey = np.sum(res_grey)
    
    # 计算绿色数据残差平方和
    x_gr, y_gr = datagreen[:, 0], datagreen[:, 1]
    res_green = (x_gr**2 + (y_gr - c_beta)**2 - R_beta**2)**2
    sum_green = np.sum(res_green)
    
    # 返回总残差平方和
    return sum_grey + sum_green

# ---------------------- 3. 设置初始猜测与变量边界 ----------------------
# 初始猜测值（需根据实际数据范围调整）
initial_guess = [4, -3, 5]  # R_beta, c_beta, R_alpha

# 变量边界：R_beta>0，c_beta<0，R_alpha>0
bounds = [(0.1, None), (None, -0.1), (0.1, None)]

# ---------------------- 4. 执行优化 ----------------------
result = minimize(objective, initial_guess, bounds=bounds, method='L-BFGS-B')

# 输出拟合结果
R_beta_fit, c_beta_fit, R_alpha_fit = result.x
c_alpha_fit = np.sqrt(R_alpha_fit**2 - R_beta_fit**2 + c_beta_fit**2)
print("拟合结果：")
print(f"绿色圆弧：R_beta={R_beta_fit:.2f}, c_beta={c_beta_fit:.2f}")
print(f"灰色圆弧：R_alpha={R_alpha_fit:.2f}, c_alpha={c_alpha_fit:.2f}")

# ---------------------- 5. 绘图验证 ----------------------
plt.figure(figsize=(8, 6))

# 绘制原始数据
plt.scatter(datagrey[:,0], datagrey[:,1], color='grey', label='灰色数据', alpha=0.6)
plt.scatter(datagreen[:,0], datagreen[:,1], color='green', label='绿色数据', alpha=0.6)

# 绘制拟合的灰色圆弧（y<0部分）
theta_grey_fit = np.linspace(np.arccos((-c_alpha_fit)/R_alpha_fit), 2*np.pi - np.arccos((-c_alpha_fit)/R_alpha_fit), 200)
x_grey_fit = R_alpha_fit * np.sin(theta_grey_fit)
y_grey_fit = c_alpha_fit + R_alpha_fit * np.cos(theta_grey_fit)
plt.plot(x_grey_fit, y_grey_fit, color='grey', linewidth=2, label='拟合灰色圆弧')

# 绘制拟合的绿色圆弧（y<0部分）
theta_green_fit = np.linspace(0, np.pi, 200)
x_green_fit = R_beta_fit * np.sin(theta_green_fit)
y_green_fit = c_beta_fit + R_beta_fit * np.cos(theta_green_fit)
plt.plot(x_green_fit, y_green_fit, color='green', linewidth=2, label='拟合绿色圆弧')

plt.xlabel('x')
plt.ylabel('y')
plt.ylim(-10, 5)
plt.legend()
plt.grid(True)
plt.show()

• 初始猜测值：需要根据实际数据的大致范围合理设置，初始值越接近真实值，优化收敛越快、结果越可靠。
• c_alpha的符号：如果拟合后灰色圆弧不符合数据分布，可尝试将c_alpha = -np.sqrt(c_alpha_sq)，具体取决于圆心的实际位置。
• 约束校验：目标函数中加入了c_alpha_sq <0时返回无穷大的判断，确保始终满足约束条件；同时通过bounds限制了c_beta<0、R_beta和R_alpha为正。

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