JuMP求解非线性问题报/(::VariableRef,::QuadExpr)未定义错误

这个错误和你有没有在目标、约束上加@NLobjective/@NLconstraint宏没有直接关系，核心问题是JuMP的非线性宏只解析直接写在宏参数块内的表达式，宏外普通Julia作用域下对JuMP变量执行的非线性运算，不会走非线性建模流程，只会调用线性/二次表达式的运算接口——而线性/二次表达式原生不支持变量除二次式、三角函数、任意次幂这类非线性操作。
你代码里的触发逻辑非常明确：

• 你在宏块外直接调用普通Julia函数hamiltonien(x, eps)、hamiltonien(y, eps)，传入的x/y是JuMP定义的优化变量数组，函数内部写的./除法、sin/cos三角函数、平方运算全都是在普通Julia作用域执行的，根本没被JuMP的非线性解析器处理，一碰到VariableRef / QuadExpr这类非线性除法就直接抛错。
• 另外你代码里还有个隐藏逻辑bug：[sin(z[i,4]) for i in size(z[:,4])]写法完全错误，size(z[:,4])返回的是数组维度元组，循环只会取到元组里的维度值（也就是数组长度），不会遍历所有行索引，就算解决了非线性报错，这里算出来的结果维度和数值都是错的。

推荐方案：用@NLexpression预定义非线性表达式，所有非线性逻辑移入宏作用域 #

不要在宏外单独写普通函数处理优化变量，所有涉及变量的非线性计算要么直接写在目标/约束宏块里，要么用@NLexpression提前定义好可复用的非线性表达式，供后续约束、目标调用。
注意JuMP非线性接口不支持Julia原生广播（.+/./这类点运算）、数组切片批量操作，最好逐索引定义，避免解析失败，参考修改：

function bvpsolve(eps,N)
    sys = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 5))
    set_optimizer_attribute(sys,"tol",1e-8)
    set_optimizer_attribute(sys,"constr_viol_tol",1e-6)

    @variables(sys, begin
                   tf       
                   x[1:N+1 , 1:n]
                   y[1:N+1 , 1:n]
                   0. ≤ h[1:N] ≤ 10
               end)

    Δt = (tf-t0)/N

    # 逐点预定义x对应的哈密顿量非线性表达式
    @NLexpression(sys, hx[i=1:N+1],
        # 宏内部支持局部变量赋值，不用在外层算u
        u = x[i,8]/(2*eps*Cd2*x[i,6]*x[i,2]);
        x[i,5]*x[i,1]*sin(x[i,4]) + x[i,6]*(
            T(x[i,1])/x[i,3] - phi(x[i,1])*S*x[i,2]^2/(2*x[i,3])*(Cd1+Cd2*u^2) - g*sin(x[i,4])
        ) - x[i,7]*Cs(x[i,2])*T(x[i,1]) + 1/eps*x[i,8]*(
            phi(x[i,1])*S*x[i,2]*u/(2*x[i,3]) - g/x[i,2]*cos(x[i,4])
        )
    )
    # 同理定义y对应的哈密顿量
    @NLexpression(sys, hy[i=1:N+1],
        u = y[i,8]/(2*eps*Cd2*y[i,6]*y[i,2]);
        y[i,5]*y[i,1]*sin(y[i,4]) + y[i,6]*(
            T(y[i,1])/y[i,3] - phi(y[i,1])*S*y[i,2]^2/(2*y[i,3])*(Cd1+Cd2*u^2) - g*sin(y[i,4])
        ) - y[i,7]*Cs(y[i,2])*T(y[i,1]) + 1/eps*y[i,8]*(
            phi(y[i,1])*S*y[i,2]*u/(2*y[i,3]) - g/y[i,2]*cos(y[i,4])
        )
    )

    # 目标函数直接写在@NLobjective内
    @NLobjective(sys, Min, 
        sum(sum((x[i,j]-y[i,j])^2 for i in 1:N+1) for j in 1:n)/N + α*sum((h[i]-Δt)^2 for i in 1:N)
    )

    # 边界约束
    @NLconstraints(sys, begin
                       con_h0, x[1,1]   - 3480. == 0     
                       con_hf, x[N+1,1] - 9144. == 0  
                       con_v0, x[1,2]   - 151.67 == 0     
                       con_vf, x[N+1,2] - 191. == 0
                       con_m0, x[1,3]   - 69000. == 0
                       con_mf, x[N+1,3] - 68100. == 0
                       con_g0, x[1,4]   - 69000. == 0
                       con_gf, x[N+1,4] - 68100. == 0                
                   end)

    # 后续配点约束也全部用@NLconstraint逐点写，直接引用上面定义的hx、hy即可，不要在宏外调用hvfun这类普通函数处理变量
    # ...
end

如果你用到的T/phi/Cs/hvfun是自定义的复杂Julia函数，不能直接在@NL宏里调用，需要先用@operator宏把函数注册到当前JuMP模型，再在宏块内使用。

• 所有和优化变量相关的非线性操作（除法、三角函数、幂运算、自定义函数调用），必须全部放在@NLobjective/@NLconstraint/@NLexpression/@operator这几个JuMP非线性宏的内部，宏外普通代码里只能做变量定义、参数赋值这类不涉及变量非线性运算的操作。
• 不要在非线性表达式里用Julia的广播语法、数组切片、数组推导式做批量计算，JuMP的非线性解析器不支持这类语法，逐索引定义是兼容性最高的写法。
• 把你原来代码里错误的数组推导式逻辑全部替换成逐索引取值，避免数值计算错误。

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