# Define a recursive JuMP.@expression

**URL:** <https://discourse.julialang.org/t/define-a-recursive-jump-expression/129848>\
**Category:** Optimization (Mathematical)\
**Tags:** question, jump\
**Created:** [June 13, 2025, 1:26am UTC](https://discourse.julialang.org/t/define-a-recursive-jump-expression/129848 "2025-06-13T01:26:06Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [June 13, 2025, 1:26am UTC](https://discourse.julialang.org/t/define-a-recursive-jump-expression/129848/1 "2025-06-13T01:26:06Z")

</div>

Is this a good syntax that can be supported?

```julia
julia> using JuMP; model = Model();

julia> x0 = 1.3; # initial state

julia> @variable(model, u[t = 1:3], Bin); # actions

julia> @expression(model, x[t = 1:3], (t == 1 ? x0 + 1.3u[t] : x[t-1] + 1.3u[t])) # states
ERROR: UndefVarError: `x` not defined in `Main`
Suggestion: check for spelling errors or missing imports.
Stacktrace:
  [1] macro expansion
    @ K:\judepot1115\packages\MutableArithmetics\tNSBd\src\rewrite_generic.jl:67 [inlined]
  [2] macro expansion
    @ K:\judepot1115\packages\MutableArithmetics\tNSBd\src\rewrite.jl:371 [inlined]
  [3] macro expansion
    @ K:\judepot1115\packages\JuMP\RGIK3\src\macros.jl:264 [inlined]
  [4] macro expansion
    @ K:\judepot1115\packages\JuMP\RGIK3\src\macros\@expression.jl:86 [inlined]
  [5] (::var"#5#6"{Model})(t::Int64)
    @ Main K:\judepot1115\packages\JuMP\RGIK3\src\Containers\macro.jl:550
  [6] #84
    @ K:\judepot1115\packages\JuMP\RGIK3\src\Containers\container.jl:85 [inlined]
  [7] iterate
    @ .\generator.jl:48 [inlined]
  [8] collect_to!(dest::Vector{…}, itr::Base.Generator{…}, offs::Int64, st::Tuple{…})
    @ Base .\array.jl:849
  [9] collect_to_with_first!(dest::Vector{…}, v1::AffExpr, itr::Base.Generator{…}, st::Tuple{…})
    @ Base .\array.jl:827
 [10] collect(itr::Base.Generator{JuMP.Containers.VectorizedProductIterator{…}, JuMP.Containers.var"#84#85"{…}})        
    @ Base .\array.jl:801
 [11] map(f::Function, A::JuMP.Containers.VectorizedProductIterator{Tuple{Base.OneTo{Int64}}})
    @ Base .\abstractarray.jl:3399
 [12] container
    @ K:\judepot1115\packages\JuMP\RGIK3\src\Containers\container.jl:85 [inlined]
 [13] container
    @ K:\judepot1115\packages\JuMP\RGIK3\src\Containers\container.jl:71 [inlined]
 [14] container(f::Function, indices::JuMP.Containers.VectorizedProductIterator{…}, ::Type{…}, names::Vector{…})        
    @ JuMP.Containers K:\judepot1115\packages\JuMP\RGIK3\src\Containers\container.jl:75
 [15] macro expansion
    @ K:\judepot1115\packages\JuMP\RGIK3\src\macros.jl:400 [inlined]
 [16] top-level scope
    @ REPL[4]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia> # since it fails, I write the following my own code

julia> x = Vector{AffExpr}(undef, length(u))
3-element Vector{AffExpr}:
 #undef
 #undef
 #undef

julia> x[1] = x0 + 1.3u[1]
1.3 u[1] + 1.3

julia> for t in 2:length(x)
           x[t] = x[t-1] + 1.3u[t]
       end

julia> x # is this good?
3-element Vector{AffExpr}:
 1.3 u[1] + 1.3
 1.3 u[1] + 1.3 u[2] + 1.3
 1.3 u[1] + 1.3 u[2] + 1.3 u[3] + 1.3

```

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [June 13, 2025, 1:40am UTC](https://discourse.julialang.org/t/define-a-recursive-jump-expression/129848/2 "2025-06-13T01:40:38Z")

</div>

What about this

```julia
julia> using JuMP; model = Model();

julia> x0 = 1.3; # initial state

julia> co = 1.3; # coefficient of actions

julia> @variable(model, u[t = 1:3], Bin); # actions

julia> function fx(x0, co, u, t)
           if t == 1
               return x0 + u[t]co
           else
               return fx(x0, co, u, t-1) + u[t]co
           end
       end
fx (generic function with 1 method)

julia> f = t -> fx(x0, co, u, t) # ready this generating function
#5 (generic function with 1 method)

julia> @expression(model, x[t = 1:3], f(t))
3-element Vector{AffExpr}:
 1.3 u[1] + 1.3
 1.3 u[1] + 1.3 u[2] + 1.3
 1.3 u[1] + 1.3 u[2] + 1.3 u[3] + 1.3

```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [June 13, 2025, 3:00am UTC](https://discourse.julialang.org/t/define-a-recursive-jump-expression/129848/3 "2025-06-13T03:00:47Z")

</div>

Recursive expressions are not supported. I’d write your expression as:

```julia
using JuMP
model = Model()
@variable(model, u[t = 1:3], Bin)
x0 = 1.3
@expression(model, x[t = 1:3], 1.3 * u[t])
add_to_expression!(x[1], x0)
for t in 2:3
    add_to_expression!(x[t], x[t-1])
end

```

Your section option is also fine, but it is going to recursively re-evaluate the intermediate expressions.

You could also just write it explicitly:

```julia
using JuMP
model = Model()
@variable(model, u[t = 1:3], Bin)
x0 = 1.3
@expression(model, x[t = 1:3], x0 + sum(1.3 * u[i] for i in 1:t))

```

Recursive expressions are usually a sign that you could extract them into a variable:

```julia
using JuMP
model = Model()
@variable(model, u[t in 1:3], Bin)
@variable(model, x[t in 1:3])
@constraint(model, [t in 1:3], x[1] == (t == 1 ? x0 : x[t-1]) + 1.3 * u[t])

```

This allows the solver to more easily exploit the structure in the problem.

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [June 13, 2025, 3:35am UTC](https://discourse.julialang.org/t/define-a-recursive-jump-expression/129848/4 "2025-06-13T03:35:35Z")

</div>

The `add_to_expression!` looks good, I’ll go with it.

> [@odow](#):
>
> allows the solver to more easily exploit the structure

How do you know it can? In that case it will have both more columns and rows.

I did a test with my application instance

 ![image](https://global.discourse-cdn.com/julialang/original/3X/5/2/52d8be2e13491fc58137a634acece2731407d202.png)

It appears that the formulation by avoiding `@variable` (i.e. using recursive expression) is a bit faster.

The code is as follows. There are 2 `if-end` blocks in the end. The 1st `if-end` block is using `@expression` only. The 2nd `if-end` block is defining additional `@variable` and `@constraint`.

> **Code**
>
> ```julia
> import SparseArrays, Random, PGLib; RawDict = PGLib.pglib("pglib_opf_case118_ieee.m"); Random.seed!(1);
> function noise() return 1e-6rand(-1:1e-9:1) end;
> function noise(N) return 1e-6rand(-1:1e-9:1, N) end;
> function perturb!(v) return for i in eachindex(v) v[i] += noise() end end;
> function fill_Dmat!(Dmat)
> for c in 2:size(Dmat, 2)
> a = rand(0.95:7e-9:1.17, D) + rand(-0.017:7e9:0.017, D)
> Dmat[:, c] = view(Dmat, :, c-1) .* a
> end
> end; function write_to!(Dmat, Darray, t)
> for c in 1:size(Dmat, 2)
> Darray[:, t, c] = view(Dmat, :, c)
> end
> end; function rectify_dir!(bft)
> for r in 1:size(bft, 1)
> if bft[r, 1] > bft[r, 2]
> bft[r, 1], bft[r, 2] = bft[r, 2], bft[r, 1]
> elseif bft[r, 1] == bft[r, 2]
> error("self-branch")
> end
> end
> end; function assert_is_connected(bft)
> B = size(bft, 1)
> pb = Vector(1:B) # primal bs
> sn = [1] # subnet
> for ite in 1:B # only du a full iteration as a conservative choice
> progress = false
> for (i, b) in enumerate(pb) # take out branch b
> if bft[b, 1] in sn
> on = bft[b, 2] # the other node
> on ∉ sn && push!(sn, on)
> elseif bft[b, 2] in sn
> on = bft[b, 1] # the other node
> on ∉ sn && push!(sn, on)
> else
> continue
> end
> popat!(pb, i); progress = true; break # fathom branch b
> end
> if progress == false
> error("ite = $ite. All the rest branches are not connected to the subnet being investigated. Check the current subnet")
> end
> 1:maximum(bft) ⊆ sn && return # The graph is proved to be connected
> end
> error("here shouldn't be reached")
> end; function bft_2_A(bft)
> B, N = size(bft, 1), maximum(bft)
> return SparseArrays.sparse([Vector(1:B); Vector(1:B)], vec(bft), [-ones(Int, B); ones(Int, B)], B, N)
> end; function adjnode(n)
> nv = Int[]
> for ci in findall(x -> x == n, bft)
> b, c = ci.I # row(= branch), col
> on = bft[b, 3 - c]
> on ∉ nv && push!(nv, on)
> end
> sort(nv)
> end; function get_b(n, m) # use this after `rectify_dir!(bft)` had been done
> n > m && ((n, m) = (m, n))
> for r in 1:size(bft, 1)
> (bft[r, 1] == n && bft[r, 2] == m) && return r
> end
> error("no such branch exist")
> end; G, G2N, Gcapinit, Gptcost = let
> G = length(RawDict["gen"])
> G2N = [RawDict["gen"]["$g"]["gen_bus"] for g in 1:G]
> cost = [RawDict["gen"]["$g"]["cost"] for g in 1:G]
> pmax = [RawDict["gen"]["$g"]["pmax"] for g in 1:G]
> gi = findall(x -> x > 1e-5, pmax)
> G2N, cost, pmax, G = G2N[gi], cost[gi], pmax[gi], length(gi)
> cost = [.001i[1] for i in cost]
> perturb!(pmax)
> perturb!(cost)
> G, G2N, pmax, cost
> end; D, D2N, Dref = let
> D = length(RawDict["load"])
> D2N = [RawDict["load"]["$d"]["load_bus"] for d in 1:D]
> Dref = [RawDict["load"]["$d"]["pd"] for d in 1:D]
> perturb!(Dref)
> D, D2N, Dref
> end;
> Y, T, Dura, Darray = let
> Y, T = 15, 3
> Dura = [6, 3, 1]
> Darray = Array{Float64, 3}(undef, D, T, Y)
> Dmat = Matrix{Float64}(undef, D, Y) # temporary
> Random.seed!(2)
> Dmiddle = rand(1.1:7e-13:1.3, D) .* Dref
> Dlarge = rand(1.4:7e-13:1.7, D) .* Dref
> Dmat[:, 1] = Dref; fill_Dmat!(Dmat); write_to!(Dmat, Darray, 1)
> Dmat[:, 1] = Dmiddle; fill_Dmat!(Dmat); write_to!(Dmat, Darray, 2)
> Dmat[:, 1] = Dlarge; fill_Dmat!(Dmat); write_to!(Dmat, Darray, 3)
> Y, T, Dura, Darray
> end;
> S, Aprob, Avail, g_hyd, Ginvcost = let
> g_hyd = 17; Gptcost[g_hyd] = 0 # 💧 pg[g = 17] is a hydro generator
> Ginvcost = 2.0Gptcost;
> Ginvcost[g_hyd] = maximum(Gptcost);
> S = 4; Aprob = 1/S * ones(S)
> Random.seed!(1)
> Avail = sin.(rand(0.37:7e-13:1.3, S, Y))
> S, Aprob, Avail, g_hyd, Ginvcost
> end; 
> (B, N) = length.((RawDict["branch"], RawDict["bus"])) # (186, 118)
> bft = let # load a rectified bft
> bft = [RawDict["branch"]["$r"][c == 1 ? "f_bus" : "t_bus"] for r in 1:B, c in 1:2]; # raw bft
> assert_is_connected(bft);
> rectify_dir!(bft);
> bft
> end; A = bft_2_A(bft);
> import JuMP, Gurobi; GRB_ENV = Gurobi.Env();
> function bpf(b, Dvec, pgvec, pfvec) # dispatch function of power flow on branch b
> b == 9 && return -pgvec[findfirst(x -> x == 10, G2N)::Int] # b8[9] = -3
> b == 7 && return bpf(9, Dvec, pgvec, pfvec) # b8[7] = -1
> b == 94 && return bpf(93, Dvec, pgvec, pfvec) # b8[94] = -1
> b == 127 && return -bpf(126, Dvec, pgvec, pfvec) # b8[127] = -1
> b == 134 && return -pgvec[findfirst(x -> x == 87, G2N)::Int] # b8[134] = -3
> b == 176 && return -pgvec[findfirst(x -> x == 111, G2N)::Int] # b8[176] = -3
> b == 184 && return Dvec[findfirst(x -> x == 117, D2N)::Int] # b8[184] = -2
> dup_bvec = [67, 76, 99, 86, 124, 139, 142]
> b ∈ dup_bvec && return false # b8[dup_bvec] .= 0; b8[dup_bvec .- 1] .= 2
> return pfvec[findfirst(x -> x == b, F2B)::Int] # [Fallback] the corresponding JuMP's decision
> end
> n8 = let # elements being
> # [1]: a bus with load only, and whose degree is larger than 2
> # [2]: a bus with gen only, and whose degree is larger than 2
> # [3]: a bus with both power injection and load
> # [-1]: a bus with load only, and whose degree is 1, we won't add KCL constr here
> # [-3]: a bus with gen only, and whose degree is 1, we won't add KCL constr here
> # [-2]: a bare bus whose degree is 2, we won't add KCL constr here
> # [0]: a bare joint bus
> n8 = zeros(Int8, N); n8[D2N] .+= 1; n8[G2N] .+= 2
> n8[117] = -1
> n8[[10, 87, 111]] .= -3
> n8[[9, 63, 81]] .= -2
> n8
> end; b8 = let # Note this procedure is specific to case118, therefore executing once is sufficient
> # elements being
> # [2]: this line represents a 2-parallel branch, works as a normal line
> # [1]: a normal line
> # [0]: this line was banned (pf ≡ 0) due to parallel redundancy, we won't allocate pf variable here
> # [-1]: power flow of this line hinges on one of its adjacent line, we won't allocate pf variable here
> # [-2]: power flow of this line hinges on the corresponding pure load, we won't allocate pf variable here
> # [-3]: power flow of this line hinges on the corresponding pure gen, we won't allocate pf variable here
> b8 = ones(Int8, size(bft, 1))
> b8[9] = -3
> b8[7] = -1
> b8[94] = -1
> b8[127] = -1
> b8[134] = -3
> b8[176] = -3
> b8[184] = -2
> dup_bvec = [67, 76, 99, 86, 124, 139, 142]
> b8[dup_bvec] .= 0; b8[dup_bvec .- 1] .= 2
> b8
> end; F, F2B, Fra = let
> F2B = findall(s -> s ∈ 1:2, b8); F = length(F2B);
> Fra = [RawDict["branch"]["$b"]["rate_a"] for b in F2B]; # f_rate_a is a used part of b_rate_a, under Dref, if Fra *= .59, then INFEASIBLE
> Fra[Fra .> 5] /= 2 # my modification
> perturb!(Fra)
> F, F2B, Fra
> end;
> Finvcost = 80abs2.([RawDict["branch"]["$b"]["br_r"] + RawDict["branch"]["$b"]["br_x"]im for b in F2B]);
> perturb!(Finvcost);
> INVG = 5ones(G) + 1e5noise(G)
> INVF = 10ones(F) + 1e5noise(F)
> Finvcost *= 5
> Ginvcost *= 5
> 
> if false # deterministic formulation with NA constr
> det = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV));
> # 1st-stage variable
> JuMP.@variable(det, invf[f = 1:F, s = 1:S], Bin);
> JuMP.@variable(det, invg[g = 1:G, y = 1:Y, s = 1:S], Bin);
> # Capacity cap
> JuMP.@expression(det, fmax[f = 1:F, s = 1:S], Fra[f] + INVF[f]invf[f, s]);
> JuMP.@expression(det, pgmax[g = 1:G, y = 1:Y, s = 1:S], invg[g, y, s]INVG[g]);
> y = 1; for g = 1:G, s = 1:S
> JuMP.add_to_expression!(pgmax[g, y, s], Gcapinit[g])
> end; for y = 2:Y, g = 1:G, s = 1:S
> JuMP.add_to_expression!(pgmax[g, y, s], pgmax[g, y-1, s])
> end;
> # operational
> JuMP.@variable(det, pf[f = 1:F, t = 1:T, y = 1:Y, s = 1:S]);
> JuMP.@variable(det, pg[g = 1:G, t = 1:T, y = 1:Y, s = 1:S] ≥ 0);
> JuMP.@constraint(det, [g = 1:G, y = 1:Y, s = 1:S], view(pg, g, :, y, s) .≤ pgmax[g, y, s]);
> g = g_hyd; JuMP.@constraint(det, [y = 1:Y, s = 1:S], # 💧
> view(pg, g, :, y, s)'Dura ≤ Avail[s, y]pgmax[g, y, s]sum(Dura));
> for f = 1:F, s = 1:S
> v, B = view(pf, f, :, :, s), fmax[f, s]
> JuMP.@constraint(det, -B .≤ v); JuMP.@constraint(det, v .≤ B)
> end
> # invest cost is w.p.1 due to NA, generation cost is scene-wise
> JuMP.@expression(det, invf_cost, Finvcost'view(invf, :, S))
> JuMP.@expression(det, invg_cost[y = 1:Y], Ginvcost'view(invg, :, y, S))
> JuMP.@expression(det, gen_cost[s = 1:S, y = 1:Y], Gptcost'view(pg, :, :, y, s)Dura)
> JuMP.@objective(det, Min, (invf_cost + sum(invg_cost)) + (Aprob'sum(gen_cost; dims = 2))[1])
> for s = 1:S-1 # NA
> JuMP.@constraint(det, view(invf, :, s) == view(invf, :, S))
> JuMP.@constraint(det, view(invg, :, :, s) == view(invg, :, :, S))
> end
> # network power balance
> for t = 1:T, y = 1:Y, s = 1:S, nState = 0:3, n ∈ findall(ø -> ø == nState, n8)
> p_gen = (nState < 2 ? false : pg[findfirst(x -> x == n, G2N)::Int, t, y, s])
> p_demand = (nState == 1 || nState == 3 ? Darray[findfirst(x -> x == n, D2N)::Int, t, y] : false)
> bv, sgv = SparseArrays.findnz(view(A, :, n))
> f = b -> bpf(b, view(Darray, :, t, y), view(pg, :, t, y, s), view(pf, :, t, y, s))
> JuMP.@constraint(det, sgv'f.(bv) + p_gen == p_demand)
> end
> JuMP.optimize!(det); JuMP.assert_is_solved_and_feasible(det; allow_local = false);
> JuMP.objective_value(det) # ✅ 9177.799530047192
> end
> 
> if true # deterministic formulation [alternative]
> det = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV));
> # 1st-stage variable
> JuMP.@variable(det, invf[f = 1:F, s = 1:S], Bin);
> JuMP.@variable(det, invg[g = 1:G, y = 1:Y, s = 1:S], Bin);
> JuMP.@variable(det, Gcapinit[g] ≤ pgmax[g = 1:G, y = 1:Y, s = 1:S] ≤ Gcapinit[g] + INVG[g]y);
> y = 1; JuMP.@constraint(det, [g = 1:G, s = 1:S], pgmax[g, y, s] == Gcapinit[g] + INVG[g]invg[g, y, s]);
> JuMP.@constraint(det, [g = 1:G, y = 2:Y, s = 1:S], pgmax[g, y, s] == pgmax[g, y-1, s] + INVG[g]invg[g, y, s]);
> JuMP.@expression(det, fmax[f = 1:F, s = 1:S], Fra[f] + INVF[f]invf[f, s]);
> # operational
> JuMP.@variable(det, pf[f = 1:F, t = 1:T, y = 1:Y, s = 1:S]); for f = 1:F, s = 1:S
> v, B = view(pf, f, :, :, s), fmax[f, s]
> JuMP.@constraint(det, -B .≤ v)
> JuMP.@constraint(det, v .≤ B)
> end
> JuMP.@variable(det, pg[g = 1:G, t = 1:T, y = 1:Y, s = 1:S] ≥ 0);
> JuMP.@constraint(det, [g = 1:G, y = 1:Y, s = 1:S], view(pg, g, :, y, s) .≤ pgmax[g, y, s]);
> g = g_hyd; JuMP.@constraint(det, [y = 1:Y, s = 1:S], # 💧
> view(pg, g, :, y, s)'Dura ≤ Avail[s, y]pgmax[g, y, s]sum(Dura));
> # invest cost is w.p.1 due to NA, generation cost is scene-wise
> JuMP.@expression(det, invf_cost, Finvcost'view(invf, :, S));
> JuMP.@expression(det, invg_cost[y = 1:Y], Ginvcost'view(invg, :, y, S));
> JuMP.@expression(det, gen_cost[s = 1:S, y = 1:Y], Gptcost'view(pg, :, :, y, s)Dura);
> JuMP.@objective(det, Min, (invf_cost + sum(invg_cost)) + (Aprob'sum(gen_cost; dims = 2))[1]);
> for s = 1:S-1 # NA
> JuMP.@constraint(det, view(invf, :, s) == view(invf, :, S))
> JuMP.@constraint(det, view(invg, :, :, s) == view(invg, :, :, S))
> end
> # real-time network power balance
> for t = 1:T, y = 1:Y, s = 1:S, nState = 0:3, n ∈ findall(ø -> ø == nState, n8)
> p_gen = (nState < 2 ? false : pg[findfirst(x -> x == n, G2N)::Int, t, y, s])
> p_demand = (nState == 1 || nState == 3 ? Darray[findfirst(x -> x == n, D2N)::Int, t, y] : false)
> bv, sgv = SparseArrays.findnz(view(A, :, n))
> f = b -> bpf(b, view(Darray, :, t, y), view(pg, :, t, y, s), view(pf, :, t, y, s))
> JuMP.@constraint(det, sgv'f.(bv) + p_gen == p_demand)
> end
> JuMP.optimize!(det); JuMP.assert_is_solved_and_feasible(det; allow_local = false);
> JuMP.objective_value(det) # ✅ 9177.799530047192
> end
> 
> ```

I guess that by using `@expression` only, we only end up with a dense “A” matrix, but the size of whom is smaller. If we use `@variable` and `@constraint`, we will end up with large and sparse “A” matrix.