方程求根局部发散问题：两相流delta_l与VSG参数稳定求解求助

需求说明 #

需要在VSL取值区间[0.001, 1]内求解目标方程得到delta_l值，再代入后续公式计算VSG，复现参考曲线图。

现有方案问题 #

• 基于scipy.optimize.fsolve、scipy.root的所有求解方法：仅当给定的delta_l初值接近真实解时，能得到VSL在[0.001, 0.1]区间的解，更换初值后全区间无有效解。
• 基于遗传算法的进化类求解方法：仅当delta_l求解范围设为[0.001, 0.065]时能得到预期结果，放宽或收窄求解范围所得结果不符合要求。
• 预期目标：将delta_l的求解范围设为[0.0001, 0.2]时，算法仍能收敛到正确解。

现有实现代码 #

import pygad
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
import math

D= 0.06  ### Diameter (m)
ROL= 997.9 #### Density liquid  Kg/m3
ROG= 1.2  ##### density gas    kg/m3
UL= 1.1*(10**(-3))### viscosit Pa.S
UG= 0.000018  ### viscosit Pa.S
G= 9.81   # gravitational acceleration  m/s2
sigma= 0.06    #### Interfacial tension  N/m
SL=3.14*D
### constants


list_of_derivatives=[]
list_of_interfacial_tension=[]
list_of_interfacial_tension_gas=[]
    
list_of_VSG=[]   ###### main values we are looking for



alpha= 20   #### inclination from horizontal   Degrees
list_of_liquid_velocities_VSL= np.linspace(0.001,1,80) #va


##################### barneaaaaaa   ##############

for VSL in list_of_liquid_velocities_VSL: 
    AL=3.14*(D**2)*(0.01-0.01**2)
    
    DL=4*AL/SL
    A=3.14*(0.5*D)**2
    VL=VSL*A/AL
    if (D*ROL*VSL/(UL))<2300:
        n=1
        CL=16
        m=0.2
        CG=0.046
        
       
    else:
        n=0.2
        CL=0.046
        m=0.2
        CG=0.046
        
        
    desired_output = 0
    
    def fitness_func(x,x_idx):
        derivative = ((G*(ROL-ROG)*D*math.sin(alpha*3.14/180)*((1-2*x)**2-2*(x-x**2)))-((CL/16)*ROL*((D*ROL/UL)**(-n))*(VSL**(2-n))*(((x-x**2)+(1-2*x)**2)/(x-x**2)**3)))

        fitness = (1/np.abs(derivative - desired_output))
        
        return fitness

    fitness_function = fitness_func

    num_generations = 100
    num_parents_mating = 4

    sol_per_pop = 8
    num_genes = 1
    init_range_low = 0.001
    init_range_high = 0.101

    parent_selection_type = "sss"
    keep_parents = 1

    crossover_type = "single_point"

    mutation_type = "random"
    mutation_percent_genes = 20

    gene_space = [{'low': 0.001, 'high': 0.065}]
    ga_instance = pygad.GA(num_generations=num_generations,
                       num_parents_mating=num_parents_mating,
                       fitness_func=fitness_func,
                       sol_per_pop=sol_per_pop,
                       num_genes=num_genes,
                       init_range_low=init_range_low,
                       init_range_high=init_range_high,
                       parent_selection_type=parent_selection_type,
                       keep_parents=keep_parents,
                       crossover_type=crossover_type,
                       mutation_type=mutation_type,
                       mutation_percent_genes=mutation_percent_genes,
                       gene_space= gene_space)


    ga_instance.run()

    x, solution_fitness, solution_idx = ga_instance.best_solution()
    
   

    derivative = ((G*(ROL-ROG)*D*np.sin(alpha*3.14/180)*((1-2*x)**2-2*(x-x**2)))-((CL/16)*ROL*((D*ROL/UL)**(-n))*(VSL**(2-n))*(((x-x**2)+(1-2*x)**2)/(x-x**2)**3))) ### should be equal to 0


    list_of_derivatives.append(abs(derivative))
    
    interfacial_tenstion=((G*(ROL-ROG)*D*np.sin(alpha*3.14/180)*(x-x**2)*(1-2*x))+((CL/32)*ROL*(D*ROL/UL)**(-n)*VSL**(2-n)*((1-2*x)/(x-x**2)**2)))  ### should be equal to interfacial_tension_gas
    VSG  =  (interfacial_tenstion*((1-2*x)**4) /((0.5)*(CG)*((D*ROG/UG)**(-m))*(1+300*x)*(ROG)))**(1/(2-m))
    
    list_of_interfacial_tension.append(interfacial_tenstion)
    
    interfacial_tension_gas= ((((VSG)**(2-m))*(0.5)*(CG)*((D*ROG/UG)**(-m))*(1+300*x)*(ROG))/((1-2*x)**4))/(G*(ROL-ROG)*D)  ### this should equal to interfacial_tension
    list_of_interfacial_tension_gas.append(interfacial_tension_gas)
    
    list_of_VSG.append(VSG)
    
    
VSLexperiment= [0.02,0.06,0.1]
VSGexperiment=[13.94,14.54,19.09]

plt.plot()
plt.xscale("log")
plt.yscale("log")
plt.xlim(1,100)
plt.ylim(0.001,1)
plt.xlabel("VSG", fontsize=9)
plt.ylabel("VSL",fontsize=9)
plt.grid(True, which="both", ls=":", color='0.002', linewidth=0.4)
plt.scatter(VSGexperiment,VSLexperiment, color="Black", linewidth="0.05", label="VsgExperiment")
plt.plot (list_of_VSG,list_of_liquid_velocities_VSL, label='VBarnea', color='blue')  

可行解决思路 #

1. 优先使用二分法求解
   这个问题属于单变量求根场景，可先对每个VSL值，在[0.0001, 0.2]区间内做100个点的粗粒度采样，计算目标函数值定位过零点的区间，只要区间两端函数符号相反，二分法可以100%收敛到根，完全不依赖初值，稳定性和效率都远高于遗传算法和梯度类求解器。
2. 优化梯度类求解器的启动逻辑
   如果要使用scipy.root类的梯度求解器，不要使用固定初值，而是采用暖启动策略：把前一个VSL点求解得到的delta_l作为下一个VSL点的初值，因为VSL是连续变化的，delta_l的取值也是连续渐变的，该策略能大幅提升全区间收敛成功率。
3. 遗传算法参数优化
   如果要保留遗传算法实现，做以下修改即可支持更大的求解范围：

• 种群规模提升到20以上，迭代次数增加到300代，避免小种群过早收敛到局部最优
• 把适应度函数从1/|f(x)|修改为-np.abs(derivative)，避免f(x)接近0时适应度值爆炸导致的收敛不稳定
• 精英保留数量提升到3，保证每一代的最优解不会被交叉、变异操作破坏

4. 增加边界惩罚规避奇点
   目标函数中存在(x-x**2)**3的分母项，x接近0或1时函数值会异常突变，可在搜索逻辑中增加边界惩罚，只要x超出物理合理范围就赋予极低的适应度，避免算法搜索到无意义的奇点区域。

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