OptimizationProblem(sys; grad = true) fails when the cost contains a lazy array reduction such as sum(A * x) nested inside another scalar operation (2 * sum(A * x), sum(abs2, D*u - p) + sum(abs2, u), etc.). The root cause is in Symbolics' gradient/expand_derivatives. It reproduces without MTK, but it affects the MTK scalar cost gradient (calculate_cost_gradient in lib/ModelingToolkitBase/src/systems/codegen.jl, which calls Symbolics.gradient(cost(sys), flat_unknowns(sys))).

I found this while working on #5196 (multi-objective optimization tests).

Minimal reproducer

using ModelingToolkitBase, Symbolics

@variables x[1:2]
A = [1.0 2.0; 3.0 4.0]

# Symbolics only
Symbolics.gradient(sum(A * x), collect(x))          # works: [4.0, 6.0]
try
    Symbolics.gradient(2 * sum(A * x), collect(x))  # errors
catch e
    println("Symbolics: ", sprint(showerror, e))
end

# ModelingToolkitBase
sys = complete(System(Equation[], [x], []; costs = [2 * sum(A * x)], name = :sys))
try
    OptimizationProblem(sys, [x => [0.1, 0.2]]; grad = true)
catch e
    println("MTK: ", sprint(showerror, e))
end

Output:

Symbolics: Differentiation of `x[1]` inside lazy array expressions is not supported; scalarize the expression with `Symbolics.scalarize` first
MTK: Differentiation of `x[1]` inside lazy array expressions is not supported; scalarize the expression with `Symbolics.scalarize` first

Symbolics.gradient(Symbolics.scalarize(expr), collect(x)) works. The multiobjective = true, jac = true path also works because Symbolics.jacobian scalarizes its inputs.

Root cause (Symbolics src/diff.jl)

• executediff has a special case for a top-level mapreduce ArrayOp (BSImpl.ArrayOp(; output_idx) && if isempty(output_idx) → SymbolicUtils.scalarize(arg)). That is why sum(A*x) alone works.
• When the reduction is an argument of +/*, executediff first calls occursin_info(D.x, arg) on the whole expression. _occursin_info then recurses into the mapreduce ArrayOp with fail = false, reaches the array-shaped A*x, and hits the error("Differentiation of ... inside lazy array expressions is not supported ...") branch before it ever reaches the scalarizing special case.
• The any(predicate, arguments(expr)) loop short-circuits, so whether the error appears depends on the hash-based order of the Add/Mul arguments. If an earlier argument already contains x[1], the lazy term is never inspected. That makes the failure look value-dependent. With D = [2.0 -1 0; -1 2 -1; 0 -1 2], sum(abs2, D*u) + sum(abs2, u) differentiates fine because mapreduce(abs2, add_sum, u) sorts first. With D = 16*[...] the D*u term sorts first and the same expression errors. The original trigger was
  System(Equation[], [u], [src]; costs = [sum(abs2, D*u - src), sum(abs2, u)]) with u[1:3] and D = 16*[2.0 -1 0; -1 2 -1; 0 -1 2].

Possible fixes: (a) in Symbolics, scalarize the scalar-valued mapreduce ArrayOp sub-terms, or have occursin_info treat them conservatively instead of erroring. (b) In MTK, as a workaround, scalarize the cost before Symbolics.gradient / Symbolics.hessian in calculate_cost_gradient (and the Hessian counterpart). The Symbolics fix is the real one. I'm filing this here because it is the user-visible failure in MTK.

Versions

ModelingToolkitBase 1.77.0 (monorepo master e167d97), Symbolics 7.41.0, SymbolicUtils 4.47.2.

Julia Version 1.11.9
Commit 53a02c0720c (2026-02-06 00:27 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 128 × AMD EPYC 7502 32-Core Processor
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, znver2)
Threads: 1 default, 0 interactive, 1 GC (on 128 virtual cores)

🤖 Filed by an AI agent (Claude Code 2.1.280, model claude-opus-5-5); local session on amdci2.julia.csail.mit.edu.