向量ODE系统参数估计优化失败，求代码修正方案

我尝试基于合成数据进行参数估计：先通过给定参数求解ODE生成数据，再添加随机噪声。使用build_loss_objective构建代价函数后传入优化器，通过result.minimizer获取估计参数。

我已实现单变量ODE的参数估计示例，能成功得到参数k和β的估计值：

using DifferentialEquations, Random, Plots
using DiffEqParamEstim, Optim

# ODE Model
function fun(G,p,t)
 k, β  = p
 #ODE
 dG_dt = k - β*G
 return dG_dt
end

# Parameter values
p = [0.8, 0.2] #"true" Parameter values
G_ini = 0.0
tspan = (0.0,10.0)
# ODE Problem
prob = ODEProblem(fun,G_ini,tspan,p)
# solving ODE Problem
sol = solve(prob,Tsit5())

# Synthetic Data
dataset = [(t,sol(t)+0.2randn()) for t in 0:0.01:10]

# Cost function and Optimization
dataTime = [d[1] for d in dataset] # time points
dataValues = [d[2] for d in dataset] # Data

cost_function = build_loss_objective(prob,Tsit5(),L2Loss(dataTime,dataValues),maxiters=10000,verbose=false)

initialGuess = ones(2) # guessed parameter values
result = optimize(cost_function, initialGuess, BFGS())

但扩展到向量ODE系统时，运行optimize函数提示失败，无法得到估计参数。相关代码如下，请问需修改哪些部分才能让代码正常运行？

using DifferentialEquations,DiffEqParamEstim,RecursiveArrayTools
using Optim,Plots

# ODE Model
n=10
function f(du,u,p,t)
 du[1,1] = dc_n = 4*p[1]*u[1,1] -u[1,1]*u[1,2]  
 du[1,2] = dc_d = p[2]*u[1,2]^2 + u[1,1]*u[1,2]

 du[n,1] = dc_n = 3*p[1]*u[n,1] -u[n,1]*u[n,2]                     
 du[n,2] = dc_d = -2*p[2]*u[n,2]^2 + u[n,1]*u[n,2]


 for i in 2:(n-1)
    du[i,1] = dc_n = p[1]*u[i-1,1] -u[i,1]*u[i+1,2]               
    du[i,2] = dc_d = -p[2]*u[i-1,2]^2 + u[i,1]*u[i+1,2]
 end

end

u0 = hcat(ones(n), ones(n))
tspan = (0.0,10.0)
p = [1.5,2.0,1.2]

prob = ODEProblem(f,u0,tspan,p)
sol = solve(prob,Tsit5())

# Synthetic Data
dataset = [(t,sol(t) .+0.2randn()) for t in 0:0.01:10]

dataTime = [d[1] for d in dataset] # 提取时间点
dataValues = [d[2] for d in dataset] # 提取实验数据
# 到这里都能正常运行

cost_function =  build_loss_objective(prob,Tsit5(),L2Loss(dataTime,dataValues),maxiters=10000,verbose=false)
initialGuess = ones(3) # 所有参数初始值设为1
result = optimize(cost_function, initialGuess, BFGS())

你的向量ODE参数估计代码存在三个核心问题，逐一修正后即可正常运行：

1. 移除未使用的冗余参数 #

真实参数p = [1.5,2.0,1.2]中的第三个参数p[3]从未在ODE函数f中被调用，这会导致参数不可识别——优化器无法通过数据估计一个完全不影响模型输出的参数，直接引发优化失败。

修正：
把真实参数和初始猜测都改为仅包含两个有效参数：

p = [1.5,2.0]
initialGuess = ones(2)

2. 统一状态变量的数组结构 #

DifferentialEquations.jl对矩阵形式的状态变量支持有限，推荐将二维矩阵转为一维向量，避免维度不匹配问题。

修正方案（转一维向量）：

• 修改初始条件为一维向量：

u0 = vcat(ones(n), ones(n))  # 转为长度2n的一维向量

• 对应修改ODE函数，将u和du视为一维向量，重新索引：

function f(du,u,p,t)
    # 第1组（n个元素）对应原u[:,1]，第2组对应原u[:,2]
    du[1] = 4*p[1]*u[1] - u[1]*u[n+1]  
    du[n+1] = p[2]*u[n+1]^2 + u[1]*u[n+1]

    du[n] = 3*p[1]*u[n] - u[n]*u[2n]                     
    du[2n] = -2*p[2]*u[2n]^2 + u[n]*u[2n]

    for i in 2:(n-1)
        du[i] = p[1]*u[i-1] - u[i]*u[n+i+1]               
        du[n+i] = -p[2]*u[n+i-1]^2 + u[i]*u[n+i+1]
    end
end

3. 噪声添加与数据维度匹配 #

生成合成数据时，randn()默认生成标量，直接用sol(t) .+0.2randn()会导致广播错误（矩阵加标量噪声），需要生成和sol(t)维度一致的噪声矩阵。

修正：

dataset = [(t,sol(t) .+ 0.2randn(size(sol(t)))) for t in 0:0.01:10]

完整修正后的代码 #

using DifferentialEquations,DiffEqParamEstim
using Optim,Plots

# ODE Model
n=10
function f(du,u,p,t)
    du[1] = 4*p[1]*u[1] - u[1]*u[n+1]  
    du[n+1] = p[2]*u[n+1]^2 + u[1]*u[n+1]

    du[n] = 3*p[1]*u[n] - u[n]*u[2n]                     
    du[2n] = -2*p[2]*u[2n]^2 + u[n]*u[2n]

    for i in 2:(n-1)
        du[i] = p[1]*u[i-1] - u[i]*u[n+i+1]               
        du[n+i] = -p[2]*u[n+i-1]^2 + u[i]*u[n+i+1]
    end
end

u0 = vcat(ones(n), ones(n))
tspan = (0.0,10.0)
p = [1.5,2.0]  # 移除未使用的第三个参数

prob = ODEProblem(f,u0,tspan,p)
sol = solve(prob,Tsit5())

# Synthetic Data
dataset = [(t,sol(t) .+ 0.2randn(size(sol(t)))) for t in 0:0.01:10]

dataTime = [d[1] for d in dataset]
dataValues = [d[2] for d in dataset]

cost_function = build_loss_objective(prob,Tsit5(),L2Loss(dataTime,dataValues),maxiters=10000,verbose=false)
initialGuess = ones(2)
result = optimize(cost_function, initialGuess, BFGS())

# 查看估计结果
println("真实参数: ", p)
println("估计参数: ", result.minimizer)

额外优化建议 #

• 如果坚持使用矩阵状态变量，可改用ODEProblem(f, ArrayPartition(u0), tspan, p)，同时确保ODE函数中正确处理ArrayPartition的结构；
• 对于高维系统，BFGS优化器可能需要调整初始步长或改用L-BFGS等更适合高维问题的优化算法；
• 可添加verbose=true查看代价函数构建和优化过程的日志，方便排查问题。

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