使用scipy least_squares约束变量避免sqrt返回虚数及参数范围

一、解决根号内负数问题（避免虚数报错） #

你之前的残差函数存在两个关键问题：一是循环判断负数时直接返回布尔值，不符合残差需为数组的要求；二是惩罚项设计不合理，无法有效引导优化器避开非法参数组合。以下是两种可行方案：

方案1：残差中添加惩罚项（推荐） #

先通过np.maximum将根号内参数限制为极小正数，避免计算平方根时出错；再对根号内为负的位置添加大惩罚，让优化器主动避开这类参数组合：

import numpy as np
from scipy.optimize import least_squares

def base(r, c, K):
    arg = 1 - (K+1)*(c*r)**2
    # 替换负数为极小正数，防止sqrt报错
    safe_arg = np.maximum(arg, 1e-8)
    return c * r**2 / (1 + np.sqrt(safe_arg))

def base_residual(toFit, r, trueSag):
    c, K = toFit
    arg = 1 - (K+1)*(c*r)**2
    
    # 计算基础残差
    residual = base(r, c, K) - trueSag
    # 给根号内负数的位置添加惩罚
    penalty = np.where(arg < 0, 1e6 * np.abs(arg), 0)
    residual += penalty
    
    return residual

方案2：变量替换转化约束 #

将根号内的约束条件转化为参数的范围约束：观察1 - (K+1)*(c*r)^2 > 0，可令(K+1)*c² = t²（t为实数），此时约束变为(t*r)^2 < 1，结合数据的r范围可进一步限制t的取值。不过这种方法会增加参数维度，实现复杂度更高，不如惩罚项方案简洁。

二、参数范围约束（-90<K<90，c>0.1或c<-0.1） #

由于c的约束是不连续区间，无法直接用least_squares的bounds参数，可通过以下两种方式处理：

方案1：残差中添加参数约束惩罚 #

对违反参数范围的情况添加大惩罚，引导优化器选择合法参数：

def base_residual(toFit, r, trueSag):
    c, K = toFit
    param_penalty = 0
    
    # K的范围约束惩罚
    if K <= -90 or K >= 90:
        param_penalty += 1e6 * np.abs(K - np.clip(K, -90, 90))
    # c的范围约束惩罚（避开-0.1到0.1区间）
    if -0.1 <= c <= 0.1:
        # 惩罚值取到最近合法区间的距离乘以系数
        param_penalty += 1e6 * min(np.abs(c - 0.1), np.abs(c + 0.1))
    
    arg = 1 - (K+1)*(c*r)**2
    safe_arg = np.maximum(arg, 1e-8)
    residual = c * r**2 / (1 + np.sqrt(safe_arg)) - trueSag
    
    # 根号内负数惩罚
    penalty = np.where(arg < 0, 1e6 * np.abs(arg), 0)
    residual += penalty
    # 将参数惩罚扩展为与残差同维度的数组
    residual += param_penalty * np.ones_like(residual)
    
    return residual

方案2：拆分两次拟合（更可靠） #

分别针对c>0.1和c<-0.1两个区间进行拟合，最终选择残差最小的结果：

# 拟合c>0.1的情况
bounds_pos_c = ([0.1, -90], [np.inf, 90])
fit_pos = least_squares(base_residual, [0.2, 0], args=(r_points, truesurf), bounds=bounds_pos_c)

# 拟合c<-0.1的情况
bounds_neg_c = ([-np.inf, -90], [-0.1, 90])
fit_neg = least_squares(base_residual, [-0.2, 0], args=(r_points, truesurf), bounds=bounds_neg_c)

# 选择最优结果
best_fit = fit_pos if fit_pos.cost < fit_neg.cost else fit_neg

三、降低初始猜测的依赖性 #

通过生成多个随机初始点，选择最优拟合结果：

num_trials = 5
best_cost = float('inf')
best_result = None

for _ in range(num_trials):
    # 随机生成符合约束的初始猜测
    K_guess = np.random.uniform(-90, 90)
    c_guess = np.random.uniform(0.1, 10) if np.random.rand() > 0.5 else np.random.uniform(-10, -0.1)
    init_guess = [c_guess, K_guess]
    
    current_fit = least_squares(base_residual, init_guess, args=(r_points, truesurf), verbose=0)
    if current_fit.cost < best_cost:
        best_cost = current_fit.cost
        best_result = current_fit

print("最优拟合参数：", best_result.x)

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