ODE参数估计中自动前向微分OptimizationFunction编写故障排查
2026-8-3
问题
使用Julia进行化学反应动力学参数的多批次估计，单批次优化成功，但扩展到两个批次时用Optimization.jl启动优化报错。代码及错误如下：
using DifferentialEquations, Plots, DiffEqParamEstim using Optimization, ForwardDiff, OptimizationOptimJL, OptimizationNLopt using Ipopt, OptimizationGCMAES, Optimisers using Random #Experimental data, species B is NOT observed in the data times = [0.0, 0.071875, 0.143750, 0.215625, 0.287500, 0.359375, 0.431250, 0.503125, 0.575000, 0.646875, 0.718750, 0.790625, 0.862500, 0.934375, 1.006250, 1.078125, 1.150000] A_obs = [1.0, 0.552208, 0.300598, 0.196879, 0.101175, 0.065684, 0.045096, 0.028880, 0.018433, 0.011509, 0.006215, 0.004278, 0.002698, 0.001944, 0.001116, 0.000732, 0.000426] C_obs = [0.0, 0.187768, 0.262406, 0.350412, 0.325110, 0.367181, 0.348264, 0.325085, 0.355673, 0.361805, 0.363117, 0.327266, 0.330211, 0.385798, 0.358132, 0.380497, 0.383051] P_obs = [0.0, 0.117684, 0.175074, 0.236679, 0.234442, 0.270303, 0.272637, 0.274075, 0.278981, 0.297151, 0.297797, 0.298722, 0.326645, 0.303198, 0.277822, 0.284194, 0.301471] #Create additional data sets for a multi data set optimization #Simple noise added to data for testing times_2 = times[2:end] .+ rand(range(-0.05,0.05,100)) P_obs_2 = P_obs[2:end] .+ rand(range(-0.05,0.05,100)) A_obs_2 = A_obs[2:end].+ rand(range(-0.05,0.05,100)) C_obs_2 = C_obs[2:end].+ rand(range(-0.05,0.05,100)) #ki = [2.78E+00, 1.00E-09, 1.97E-01, 3.04E+00, 2.15E+00, 5.27E-01] #Target optimized parameters ki = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1] #Initial guess of parameters IC = [1.0, 0.0, 0.0, 0.0] #Initial condition for each species tspan1 = (minimum(times),maximum(times)) #tuple timespan of data set 1 tspan2 = (minimum(times_2),maximum(times_2)) #tuple timespan of data set 2 # data = VectorOfArray([A_obs,C_obs,P_obs])' data = vcat(A_obs',C_obs',P_obs') #Make multidimensional array containing all observed data for dataset1, transpose to match shape of ODEProblem output data2 = vcat(A_obs_2',C_obs_2',P_obs_2') #Make multidimensional array containing all observed data for dataset2, transpose to match shape of ODEProblem output #make dictionary containing data, time, and initial conditions keys1 = ["A","B"] keys2 = ["time","obs","IC"] entryA =[times,data,IC] entryB = [times_2, data2,IC] nest=[Dict(zip(keys2,entryA)),Dict(zip(keys2,entryB))] exp_dict = Dict(zip(keys1,nest)) #data dictionary #rate equations in power law form r = k [A][B] function rxn(x, k) A = x[1] B = x[2] C = x[3] P = x[4] k1 = k[1] k2 = k[2] k3 = k[3] k4 = k[4] k5 = k[5] k6 = k[6] r1 = k1 * A r2 = k2 * A * B r3 = k3 * C * B r4 = k4 * A r5 = k5 * A r6 = k6 * A * B return [r1, r2, r3, r4, r5, r6] #returns reaction rate of each equation end #Mass balance differential equations function mass_balances(di,x,args,t) k = args r = rxn(x, k) di[1] = - r[1] - r[2] - r[4] - r[5] - r[6] #Species A di[2] = + r[1] - r[2] - r[3] - r[6] #Species B di[3] = + r[2] - r[3] + r[4] #Species C di[4] = + r[3] + r[5] + r[6] #Species P end function ODESols(time,uo,parms) time_init = (minimum(time),maximum(time)) prob = ODEProblem(mass_balances,uo,time_init,parms) sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8,save_idxs = [1,3,4],saveat=time) #Integrate prob return sol end function cost_function(data_dict,parms) res_dict = Dict(zip(keys(data_dict),[0.0,0.0])) for key in keys(data_dict) pred = ODESols(data_dict[key]["time"],data_dict[key]["IC"],parms) loss = L2Loss(data_dict[key]["time"],data_dict[key]["obs"]) err = loss(pred) res_dict[key] = err end residual = sum(res_dict[key] for key in keys(res_dict)) @show typeof(residual) return residual end lb = [0.0,0.0,0.0,0.0,0.0,0.0] #parameter lower bounds ub = [10.0,10.0,10.0,10.0,10.0,10.0] #parameter upper bounds optfun = Optimization.OptimizationFunction(cost_function,Optimization.AutoForwardDiff()) optprob = Optimization.OptimizationProblem(optfun,exp_dict, ki,lb=lb,ub=ub,reltol=1E-8) #Set up optimization problem optsol=solve(optprob, BFGS(),maxiters=10000) #Solve optimization problem println(optsol.u) #print solution
错误信息：
ERROR: MethodError: no method matching ForwardDiff.GradientConfig(::Optimization.var"#89#106"{OptimizationFunction{true, Optimization.AutoForwardDiff{nothing}, typeof(cost_function), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float64}}, ::Dict{String, Dict{String, Array{Float64}}}, ::ForwardDiff.Chunk{2})
排查后怀疑是cost_function对ForwardDiff兼容性问题，但无法定位具体点，需要协助解决。
解决方案
核心错误修正
Optimization.jl要求代价函数的第一个参数必须是待优化的参数，第二个参数才是静态数据。原代码中cost_function(data_dict, parms)的参数顺序完全颠倒，导致ForwardDiff无法对参数求导，这是报错的直接原因。
其他细节修正
- 生成测试数据时，
rand(range(-0.05,0.05,100))的长度要与目标数组一致，否则会出现维度不匹配错误，改为rand(length(times[2:end])) .* 0.1 .- 0.05更简洁可靠。 - 简化实验数据的存储结构，用数组替代嵌套字典，减少不必要的复杂度。
修改后的完整代码
using DifferentialEquations, Plots, DiffEqParamEstim using Optimization, ForwardDiff, OptimizationOptimJL, OptimizationNLopt using Ipopt, OptimizationGCMAES, Optimisers using Random # 实验数据，物种B未观测 times = [0.0, 0.071875, 0.143750, 0.215625, 0.287500, 0.359375, 0.431250, 0.503125, 0.575000, 0.646875, 0.718750, 0.790625, 0.862500, 0.934375, 1.006250, 1.078125, 1.150000] A_obs = [1.0, 0.552208, 0.300598, 0.196879, 0.101175, 0.065684, 0.045096, 0.028880, 0.018433, 0.011509, 0.006215, 0.004278, 0.002698, 0.001944, 0.001116, 0.000732, 0.000426] C_obs = [0.0, 0.187768, 0.262406, 0.350412, 0.325110, 0.367181, 0.348264, 0.325085, 0.355673, 0.361805, 0.363117, 0.327266, 0.330211, 0.385798, 0.358132, 0.380497, 0.383051] P_obs = [0.0, 0.117684, 0.175074, 0.236679, 0.234442, 0.270303, 0.272637, 0.274075, 0.278981, 0.297151, 0.297797, 0.298722, 0.326645, 0.303198, 0.277822, 0.284194, 0.301471] # 生成多批次测试数据，添加噪声 n_points = length(times[2:end]) times_2 = times[2:end] .+ rand(n_points) .* 0.1 .- 0.05 # 生成-0.05到0.05的噪声 P_obs_2 = P_obs[2:end] .+ rand(n_points) .* 0.1 .- 0.05 A_obs_2 = A_obs[2:end] .+ rand(n_points) .* 0.1 .- 0.05 C_obs_2 = C_obs[2:end] .+ rand(n_points) .* 0.1 .- 0.05 # 参数初始值与边界 ki = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1] # 初始猜测 lb = zeros(6) # 参数下界 ub = fill(10.0, 6) # 参数上界 IC = [1.0, 0.0, 0.0, 0.0] # 初始条件 # 整理实验数据为数组结构（每个元素是一个批次的(time, obs, IC)） exp_data = [ (times, vcat(A_obs', C_obs', P_obs'), IC), (times_2, vcat(A_obs_2', C_obs_2', P_obs_2'), IC) ] # 反应速率方程 function rxn(x, k) A, B, C, P = x k1, k2, k3, k4, k5, k6 = k r1 = k1 * A r2 = k2 * A * B r3 = k3 * C * B r4 = k4 * A r5 = k5 * A r6 = k6 * A * B return [r1, r2, r3, r4, r5, r6] end # 质量平衡微分方程 function mass_balances(di, x, k, t) r = rxn(x, k) di[1] = -r[1] - r[2] - r[4] - r[5] - r[6] # A的变化率 di[2] = r[1] - r[2] - r[3] - r[6] # B的变化率 di[3] = r[2] - r[3] + r[4] # C的变化率 di[4] = r[3] + r[5] + r[6] # P的变化率 end # 求解ODE并返回预测值 function ODESols(time, uo, parms) tspan = (minimum(time), maximum(time)) prob = ODEProblem(mass_balances, uo, tspan, parms) sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8, save_idxs=[1,3,4], saveat=time) return sol end # 代价函数：参数在前，数据在后 function cost_function(parms, exp_data) total_loss = 0.0 for (time, obs, ic) in exp_data pred = ODESols(time, ic, parms) loss = L2Loss(time, obs) total_loss += loss(pred) end return total_loss end # 设置优化问题并求解 optfun = Optimization.OptimizationFunction(cost_function, Optimization.AutoForwardDiff()) optprob = Optimization.OptimizationProblem(optfun, ki, exp_data, lb=lb, ub=ub) optsol = solve(optprob, BFGS(), maxiters=10000) println("优化后的参数：", optsol.u)
相关产品推荐
相关产品推荐
