# Efficient way to implement multi-index nonlinear constraints in ExaModels? (PCE-expanded AC power flow)
**URL:** https://discourse.julialang.org/t/efficient-way-to-implement-multi-index-nonlinear-constraints-in-examodels-pce-expanded-ac-power-flow/133934
**Category:** Optimization (Mathematical)
**Tags:** question, jump, autodiff
**Created:** [November 17, 2025, 12:51pm UTC](https://discourse.julialang.org/t/efficient-way-to-implement-multi-index-nonlinear-constraints-in-examodels-pce-expanded-ac-power-flow/133934 "2025-11-17T12:51:48Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@yanlingogo](https://discourse.julialang.org/u/yanlingogo)
#### Post date: [November 17, 2025, 12:51pm UTC](https://discourse.julialang.org/t/efficient-way-to-implement-multi-index-nonlinear-constraints-in-examodels-pce-expanded-ac-power-flow/133934/1 "2025-11-17T12:51:48Z")
Hi all,
I’m implementing a PCE-expanded AC power flow model. In JuMP, one of my nonlinear constraints looks like:
```julia-auto
@constraint(opf, pfP[i=1:n_bus, k=1:K],
p[i,k] * T2_diag[k] ==
sum(
(G[i,j]*(E[i,l1]*E[j,l2] + F[i,l1]*F[j,l2])
+ B[i,j]*(F[i,l1]*E[j,l2] - E[i,l1]*F[j,l2])) * val
for (l1, l2, val) in T3_zip[k], j in 1:n_bus
)
)
```
This works well in JuMP and converts to ExaModels quickly.
However, when I try to implement the _same_ constraint directly in ExaModels (using generators, flattened tensor data, or custom kernel functions), the model construction becomes **much slower** and the resulting nnzj/nnzh are **3–7× larger** than the JuMP-generated ExaModel.
So my questions are:
- Is there a recommended pattern for writing large multi-index nonlinear sums in ExaModels?
- Why does JuMP → ExaModels generate a much smaller expression graph than my hand-written version?
- Are there examples or documentation on expressing tensor-structured nonlinear constraints efficiently in ExaModels?
Any advice or example code would be greatly appreciated.
Thanks!
---
### Author:  [@Paul\_Schrimpf](https://discourse.julialang.org/u/Paul_Schrimpf)
#### Post date: [November 17, 2025, 2:43pm UTC](https://discourse.julialang.org/t/efficient-way-to-implement-multi-index-nonlinear-constraints-in-examodels-pce-expanded-ac-power-flow/133934/2 "2025-11-17T14:43:22Z")
With ExaModels, model construction will generally be faster with fewer `constraint` and `objective` calls. If you have something like
```julia-auto
for i in 1:ncons
constraint(m, x[i]*A[i])
end
```
model construction (and first jacobian & hessian evaluation) will be faster if it is rewritten with a single `constraint` call. There is some [discussion here](https://github.com/exanauts/ExaModels.jl/discussions/185#discussioncomment-14352953).
The `nnzh` reported by ExaModels is misleading as the Hessian is stored in a non-compact format. The Hessian is stored as a vector of values and two vectors of indices with
```julia-auto
H[i,j] = sum(hvals[k] for k in eachindex(hvals) if hrow[k]==i && hcol[k]==j)
```
`hrow` and `hcol` can and often do contain repeats. The reported `nnzh` is simply the length of `hvals` and not the true number of non-zeros in the Hessian. I think `nnzj` works similarly. The Jacobian is definitely stored in the same format, but I haven’t checked that `nnzj` is simply the length of `jvals`.
I’m not sure how best to write large multi-index sums. I started a related [discussion](https://github.com/exanauts/ExaModels.jl/discussions/122), and there are two related issues, [102](https://github.com/exanauts/ExaModels.jl/issues/102) and [167](https://github.com/exanauts/ExaModels.jl/issues/167).
You might get some more concrete advice if you share your code or, better yet, a small working example.
---
### Author:  [@yanlingogo](https://discourse.julialang.org/u/yanlingogo)
#### Post date: [November 18, 2025, 9:19am UTC](https://discourse.julialang.org/t/efficient-way-to-implement-multi-index-nonlinear-constraints-in-examodels-pce-expanded-ac-power-flow/133934/3 "2025-11-18T09:19:08Z")
Thanks a lot, I’m not familiar with ExaModel, and in its documentatio I didn’t see example that is similar to my problem. so I test the code provided by the gpt:
```julia
constraint(
core,
(
let i = i, k = k,
l1_arr = t3.l1,
l2_arr = t3.l2,
val_arr = t3.val,
p_start = p_start,
p_end = p_end
# active injection at bus i and chaos index k
p_ik = sum(Cgpar[i, g] * Pgen[g, k] for g in 1:n_gen) - Pdpar[i, k]
rhs = sum(
(
Gpar[i, j] * (E[i, l1_arr[p]] * E[j, l2_arr[p]] +
F[i, l1_arr[p]] * F[j, l2_arr[p]]) +
Bpar[i, j] * (F[i, l1_arr[p]] * E[j, l2_arr[p]] -
E[i, l1_arr[p]] * F[j, l2_arr[p]])
) * val_arr[p]
for p in p_start:p_end, j in 1:n_bus;
init = 0.0
)
p_ik * T2par[k] - rhs
end
for i in 1:n_bus
);
lcon = 0.0,
ucon = 0.0
)
```
So I compare the execution time of ExaModel(JuMPmodel) and ExaModel with gpt-provieded code, the first one takes around 10s, but the second needs around 1 minute. I will check the link you give, thank you again!
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [November 23, 2025, 7:55pm UTC](https://discourse.julialang.org/t/efficient-way-to-implement-multi-index-nonlinear-constraints-in-examodels-pce-expanded-ac-power-flow/133934/4 "2025-11-23T19:55:00Z")
Hi @yanlingogo, welcome to the forum 😄
Sorry it took me a few days to reply, I was quite busy last week with [JuMP-dev 2025](https://jump.dev/meetings/jumpdev2025/).
Your JuMP code looks good, but I’m not an expert at ExaModels, so this is a question for one of the gang of @amontoison, @sshin23, and @frapac.
---
### Author:  [@sshin23](https://discourse.julialang.org/u/sshin23)
#### Post date: [November 24, 2025, 4:24am UTC](https://discourse.julialang.org/t/efficient-way-to-implement-multi-index-nonlinear-constraints-in-examodels-pce-expanded-ac-power-flow/133934/5 "2025-11-24T04:24:44Z")
The best way to handle summations within constraint is to use `constraint!`. A tutorial can be found here:
> **[Tutorial 1: Power Flow · Powertech tutorial](https://madsuite.org/powertech2025/dev/1-powerflow/#Implementing-the-power-flow-equations-with-ExaModels)**
>
> Documentation for Powertech tutorial.
This syntax allows ExaModels.jl to break a constraint potentially with a long summation into many, simple, and embarassingly parallelizable expressions. You can first create a constraint with typical `constraint`, without summation, and then, modify the constraints by adding the terms (there, you specify the index of the constraint you’re modifying and terms you’re adding)