### Solve Optimization Problem with IPOPT

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/parameters.md

This snippet shows the initial setup and execution of the IPOPT solver with default parameters. It prints the objective value of the found solution.

```julia
set_parameter!(c_param, θ, [200.0, 1.0])  # Double the penalty coefficient
result2 = ipopt(m_param)
println("Modified penalty objective: $(result2.objective)")
```

--------------------------------

### Naive Dynamics Function Example

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/constraint_augmentation.md

This example illustrates a naive approach to defining dynamics functions where each state equation has a separate, potentially complex, expression pattern. This can lead to a high number of patterns, impacting compilation time and GPU performance.

```julia
# Naive: 20 different expression patterns (one per state equation)
function dyn_1(x, u, v, t, n)
    return -2.0 * x[11,t,n] +
           1.414 * u[1,t] * u[2,t] * x[12,t,n] * cos(v[6]) + ...
end

function dyn_2(x, u, v, t, n)
    # ... different structure
end
```

--------------------------------

### Parameters in Constraints

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Shows how parameters can be incorporated into constraints, specifically for applications like warm-starting or parameter estimation. This example defines a target value for variables.

```julia
# Parameters in constraints (warm-starting, parameter estimation)
c = ExaCore(concrete = Val(true))
@add_var(c, y, 5)
@add_par(c, target, ones(5))
@add_con(c, y[i] - target[i] for i = 1:5; lcon = -0.1, ucon = 0.1)
```

--------------------------------

### Power Flow Example Pattern

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Sets up a power flow model by defining variables for voltage magnitude, power flow, and power generation. It then adds constraints for load and shunt power, arc flows, and generation.

```julia
c = ExaCore(concrete = Val(true))
n_bus = 100
@add_var(c, vm, n_bus; start = 1.0, lvar = 0.9, uvar = 1.1)
@add_var(c, pflow, 200)
@add_var(c, pgen, 20)

bus_data = [(i = i, pd = rand(), gs = 0.01) for i in 1:n_bus]
arc_data = [(i = i, bus = mod1(i, n_bus)) for i in 1:200]
gen_data = [(i = i, bus = mod1(i*5, n_bus)) for i in 1:20]

# Base: load and shunt power
@add_con(c, p_balance, b.pd + b.gs * vm[b.i]^2 for b in bus_data)

# Augment with arc flows and generation
@add_con!(c, p_balance[a.bus] += pflow[a.i] for a in arc_data)
@add_con!(c, p_balance[g.bus] += -pgen[g.i] for g in gen_data)
```

--------------------------------

### Get and Set Variable/Constraint Bounds

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Demonstrates how to retrieve and modify variable and constraint bounds and initial values after model creation using getter and setter functions. The model is then solved with Ipopt.

```julia
using ExaModels, NLPModelsIpopt

c = ExaCore(concrete = Val(true))
@add_var(c, x, 5; start = 1.0, lvar = -10.0, uvar = 10.0)
@add_con(c, g, x[i]^2 for i = 1:5; lcon = 0.0, ucon = 5.0)
@add_obj(c, (x[i] - 3)^2 for i = 1:5)

m = ExaModel(c)

# Get current values
x_start = get_start(m, x)
x_lower = get_lvar(m, x)
x_upper = get_uvar(m, x)
g_lower = get_lcon(m, g)
g_upper = get_ucon(m, g)

println("Initial x0: $x_start")
println("Variable bounds: [$x_lower, $x_upper]")

# Modify bounds in-place
set_start!(m, x, [2.0, 2.0, 2.0, 2.0, 2.0])
set_lvar!(m, x, [-5.0, -5.0, -5.0, -5.0, -5.0])
set_uvar!(m, x, [5.0, 5.0, 5.0, 5.0, 5.0])
set_lcon!(m, g, [-1.0, -1.0, -1.0, -1.0, -1.0])
set_ucon!(m, g, [10.0, 10.0, 10.0, 10.0, 10.0])

# Solve with modified bounds
result = ipopt(m; print_level = 0)
```

--------------------------------

### Add Constraints with Parameters

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/parameters.md

Define constraints that can involve both variables and parameters. This example shows a complex constraint using variables and standard mathematical functions.

```julia
@add_constraint(
    c_param,
    3x_p[i+1]^3 + 2 * x_p[i+2] - 5 +
    sin(x_p[i+1] - x_p[i+2])sin(x_p[i+1] + x_p[i+2]) +
    4x_p[i+1] - x_p[i]exp(x_p[i] - x_p[i+1]) - 3 for i = 1:(N-2)
)
```

--------------------------------

### ExaModels v0.9 vs v0.10 Model Building

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/upgrade.md

Compares the v0.9 mutable ExaCore API with the v0.10 functional and macro styles using immutable ExaCore. Both v0.10 examples require importing ExaModels.

```julia
# ── v0.9 ────────────────────────────────────────────────────────────────────
using ExaModels

n = 100
c = ExaCore()
x = variable(c, n; lvar = -1.0, uvar = 1.0, start = 0.0)
θ = parameter(c, ones(n))
s = subexpr(c, θ[i] * x[i]^2 for i in 1:n)
g = constraint(c, x[i] + x[i+1] for i in 1:n-1; lcon = -1.0, ucon = 1.0)
constraint!(c, g, i => sin(x[i+1]) for i in 1:n-1)
objective(c, s[i] for i in 1:n)
m = ExaModel(c)
```

```julia
# ── v0.10 (functional) ───────────────────────────────────────────────────────
using ExaModels

n = 100
c = ExaCore(concrete = Val(true))
c, x = add_var(c, n; lvar = -1.0, uvar = 1.0, start = 0.0)
c, θ = add_par(c, ones(n))
c, s = add_expr(c, θ[i] * x[i]^2 for i in 1:n)
c, g = add_con(c, x[i] + x[i+1] for i in 1:n-1; lcon = -1.0, ucon = 1.0)
c, _ = add_con!(c, g, i => sin(x[i+1]) for i in 1:n-1)
c, _ = add_obj(c, s[i] for i in 1:n)
m = ExaModel(c)
```

```julia
# ── v0.10 (macro) ────────────────────────────────────────────────────────────
using ExaModels

n = 100
c = ExaCore(concrete = Val(true))
@add_var(c, x, n; lvar = -1.0, uvar = 1.0, start = 0.0)
@add_par(c, θ, ones(n))
@add_expr(c, s, θ[i] * x[i]^2 for i in 1:n)
@add_con(c, g, x[i] + x[i+1] for i in 1:n-1; lcon = -1.0, ucon = 1.0)
@add_con!(c, g, i => sin(x[i+1]) for i in 1:n-1)
@add_obj(c, s[i] for i in 1:n)
m = ExaModel(c)
```

--------------------------------

### Multi-dimensional Subexpressions

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Illustrates the creation and usage of multi-dimensional subexpressions. This example defines an 'energy' subexpression based on a 2D variable 'u'.

```julia
# Multi-dimensional subexpressions
T, K = 5, 3
c = ExaCore(concrete = Val(true))
@add_var(c, u, 1:T, 1:K)

itr = [(t, k) for t in 1:T, k in 1:K]
c, energy = add_expr(c, u[t,k]^2 + 0.5*u[t,k] for (t,k) in itr)

# Use energy[t,k] in constraints
@add_con(c, energy[t,k] for (t,k) in itr; ucon = 1.0)
```

--------------------------------

### CPU Acceleration with KernelAbstractions

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Define a model for multi-threaded CPU execution using KernelAbstractions.jl. Ensure Julia is started with multiple threads (e.g., 'julia -t 4').

```julia
using ExaModels, KernelAbstractions

# Multi-threaded CPU (start Julia with: julia -t 4)
function create_model_cpu(N)
    c = ExaCore(; backend = CPU(), concrete = Val(true))
    @add_var(c, x, N; start = (mod(i,2) == 1 ? -1.2 : 1.0 for i = 1:N))
    @add_obj(c, 100*(x[i-1]^2 - x[i])^2 + (x[i-1] - 1)^2 for i = 2:N)
    @add_con(c, 3x[i+1]^3 + 2*x[i+2] - 5 + 4x[i+1] - x[i]*exp(x[i]-x[i+1]) - 3
             for i = 1:(N-2))
    return ExaModel(c)
end

m_cpu = create_model_cpu(1000)
```

--------------------------------

### Augmenting Constraints with add_con!

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/constraint_augmentation.md

Use `add_con!` to incrementally build constraints by adding terms. Each call contributes to a single expression pattern, which is essential for performance, especially on GPUs. This example shows augmenting a base constraint with contributions from arcs and generators.

```julia
c = ExaCore(concrete = Val(true))
@add_var(c, x, n)

# Create the base constraint (e.g., a balance equation)
@add_con(c, base, rhs[i] for i in 1:n_buses)

# Augment with contributions from different sources
@add_con!(c, base, a.bus => flow[a.i] for a in arcs)
@add_con!(c, base, g.bus => -gen[g.i] for g in generators)
```

--------------------------------

### Run Local Benchmarks

Source: https://github.com/exanauts/examodels.jl/blob/main/benchmark/README.md

Use 'make' commands to run benchmarks locally. Specify backends like 'cuda', 'amdgpu', 'oneapi', 'metal', or 'full' for all available backends. 'make compare' re-runs comparisons without re-benchmarking.

```sh
# CPU only (default)
make

# Single backend
make cuda
make amdgpu
make oneapi
make metal

# All available backends
make full

# Re-run comparison without re-benchmarking
make compare
```

--------------------------------

### Create ExaModel from ExaCore

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/parameters.md

Instantiate an ExaModel using the configured ExaCore. This model can then be solved using optimization solvers.

```julia
m_param = ExaModel(c_param)
```

--------------------------------

### Create ExaModel and Solve with Ipopt

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Converts an ExaCore object into a solvable NLPModel and then solves it using the Ipopt solver. This is the standard workflow for optimizing models defined with ExaModels.jl.

```julia
using ExaModels, NLPModelsIpopt

# Build the model
c = ExaCore(concrete = Val(true))
@add_var(c, x, 10; start = (mod(i,2) == 1 ? -1.2 : 1.0 for i = 1:10))
@add_obj(c, 100*(x[i-1]^2 - x[i])^2 + (x[i-1] - 1)^2 for i = 2:10)
@add_con(c,
    3x[i+1]^3 + 2*x[i+2] - 5 + sin(x[i+1]-x[i+2])*sin(x[i+1]+x[i+2]) +
    4x[i+1] - x[i]*exp(x[i]-x[i+1]) - 3 for i = 1:8
)

# Create the ExaModel
m = ExaModel(c)

# Solve with Ipopt
result = ipopt(m; print_level = 0)
println("Status: $(result.status)")
println("Objective: $(result.objective)")
println("Iterations: $(result.iter)")
```

--------------------------------

### Power Balance Constraints with add_con!

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/constraint_augmentation.md

Decompose complex power balance constraints into simple, uniform patterns using `add_con!`. This approach results in fewer, more efficient GPU kernels compared to a naive implementation.

```julia
# Base: load and shunt (one pattern, iterated over buses)
@add_con(c, p_balance, b.pd + b.gs * vm[b.i]^2 for b in bus_data)

# Arc flows (one pattern, iterated over arcs)
@add_con!(c, p_balance[a.bus] += p[a.i] for a in arc_data)

# Generation (one pattern, iterated over generators)
@add_con!(c, p_balance[g.bus] += -pg[g.i] for g in gen_data)
```

--------------------------------

### Benchmark Output Format

Source: https://github.com/exanauts/examodels.jl/blob/main/benchmark/README.md

The 'runbenchmark.jl' script outputs a table with performance metrics. 'compare.jl' shows relative timing ratios, where values less than 1.0 indicate improvements.

```text
==========================================================================================
  name                        nvar   ncon |      obj     cons     grad      jac     hess
==========================================================================================
  rosenrock-1000-CUDA         2994   2996 | 1.23e-05 4.56e-06 2.34e-05 3.45e-06 6.78e-06
  ...

```

```text
Relative timing: current / main  (values < 1.0 are improvements)

==============================================================
  name                        |      obj     cons     grad      jac     hess
==============================================================
  rosenrock-1000-CUDA         |    0.982    1.003    0.971    0.995    0.988
  ...

```

--------------------------------

### Testing ExaModels Extension with NLPModelsIpopt

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/develop.md

This Julia script defines a test set for the LuksanVlcekModels extension. It creates a model, solves it using Ipopt, and asserts the status, solution, and multipliers against expected values.

```julia
# test/runtest.jl

using Test, LuksanVlcekModels, NLPModelsIpopt

@testset "LuksanVlcekModelsTest" begin
    m = luksan_vlcek_model(10)
    result = ipopt(m)

    @test result.status == :first_order
    @test result.solution ≈ [
        -0.9505563573613093
        0.9139008176388945
        0.9890905176644905
        0.9985592422681151
        0.9998087408802769
        0.9999745932450963
        0.9999966246997642
        0.9999995512524277
        0.999999944919307
        0.999999930070643
    ]
    @test result.multipliers ≈ [
        4.1358568305002255
        -1.876494903703342
        -0.06556333356358675
        -0.021931863018312875
        -0.0019537261317119302
        -0.00032910445671233547
        -3.8788212776372465e-5
        -7.376592164341867e-6
    ]
end
```

--------------------------------

### Add Parameters using Macro Syntax

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Illustrates adding parameters using the @add_par macro, providing a name and an array of values. This offers a more concise way to define parameters.

```julia
@add_par(c, penalty, [200.0, 50.0])
```

--------------------------------

### Solve ExaModel with Ipopt

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/parameters.md

Solve the created ExaModel using the Ipopt solver. The result includes the optimal objective value.

```julia
result1 = ipopt(m_param)
println("Original objective: $(result1.objective)")
```

--------------------------------

### Solve with MadNLP

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Demonstrates solving the ExaModel using the MadNLP solver, highlighting its capability to support GPU acceleration.

```julia
# Solve with MadNLP (also supports GPU)
using MadNLP
result = madnlp(m)
```

--------------------------------

### Create ExaCore Model Container

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Instantiate an ExaCore container with specified precision, backend, and immutability.

```julia
using ExaModels

# Create a basic ExaCore with Float64 (default)
c = ExaCore(concrete = Val(true))

# Create with Float32 precision
c = ExaCore(Float32; concrete = Val(true))

# Create with GPU backend (requires CUDA.jl)
using CUDA, KernelAbstractions
c = ExaCore(Float32; backend = CUDABackend(), concrete = Val(true))

# Create with multi-threaded CPU backend
using KernelAbstractions
c = ExaCore(; backend = CPU(), concrete = Val(true))
```

--------------------------------

### ExaModels v0.10 Core API Usage

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/upgrade.md

Demonstrates the recommended functional style for building models with the immutable ExaCore in v0.10. Requires importing ExaModels.

```julia
using ExaModels

c = ExaCore(concrete = Val(true))
c, x = add_var(c, 10; lvar = 0.0)
c, _  = add_obj(c, x[i]^2 for i in 1:10)
m = ExaModel(c)
```

--------------------------------

### Initialize ExaCore with Concrete Values

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/parameters.md

Initialize an ExaCore with concrete values to ensure type stability and performance. This is the first step in building a parametric model.

```julia
using ExaModels, NLPModelsIpopt
c_param = ExaCore(concrete = Val(true))
```

--------------------------------

### Rebinding Core in ExaModels

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/patterns.md

Demonstrates how ExaCore is immutable and each add operation returns a new core, requiring rebinding.

```julia
c = ExaCore(concrete = Val(true))
c, x = add_var(c, 10)       # c is rebound to a new ExaCore
c, _ = add_obj(c, x[i]^2 for i in 1:10)  # c is rebound again
```

--------------------------------

### Returning Core for Further Composition

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/patterns.md

Demonstrates how to return the core from a function to allow for further model building in stages across multiple functions.

```julia
function add_dynamics!(c, x, u, data)
    @add_con(c, x[i+1] - x[i] - u[i] * data[i].dt for i in 1:length(data))
    return c   # <-- return updated core for further building
end

function full_model()
    c = ExaCore(concrete = Val(true))
    @add_var(c, x, 11)
    @add_var(c, u, 10)
    c = add_dynamics!(c, x, u, data)   # rebind c with the returned core
    @add_obj(c, x[i]^2 for i in 1:11)
    return ExaModel(c)
end
```

--------------------------------

### Advanced Expression Building with Constant, exa_sum, exa_prod

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Utilize `Constant` for compile-time optimized values and `exa_sum`/`exa_prod` for programmatic construction of sum and product nodes over ranges. Useful for building complex expressions dynamically.

```julia
using ExaModels

# Constant nodes are optimized at compile time
# x * Constant(1) → x
# x + Constant(0) → x
# x ^ Constant(2) → abs2(x)

c = ExaCore(concrete = Val(true))
@add_var(c, x, 5)

# Manual sum/product over a range using exa_sum
# Useful when building expressions programmatically
using ExaModels: exa_sum, Constant

# Sum of squares: equivalent to sum(x[k]^2 for k = 1:3)
@add_obj(c, exa_sum(k -> x[k]^2, Val(1:3)) for i = 1:5)

# Product example
@add_con(c, exa_prod(k -> x[k], Val(1:3)) for i = 1:5; ucon = 1.0)
```

--------------------------------

### Add Parameters with Dimensions

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Shows how to add parameters with specific dimensions, initializing them with a default value. This is useful for parameters that have a fixed size.

```julia
c, params = add_par(c, 10; value = 1.0)
```

--------------------------------

### Extract Constraint Multipliers (Dual Values)

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Demonstrates how to obtain the dual values (Lagrange multipliers) associated with the model's constraints using the `multipliers` function.

```julia
# Get constraint multipliers (dual values)
g_mult = multipliers(result, g)
println("Constraint multipliers: $g_mult")
```

--------------------------------

### Run Manual Benchmarks

Source: https://github.com/exanauts/examodels.jl/blob/main/benchmark/README.md

Manually run benchmarks for a specific revision and backend using 'runbenchmark.jl'. Compare results using 'compare.jl'. Ensure the project is set up with 'julia --project=.'

```sh
# Build results for a specific revision + backend
julia runbenchmark.jl main    cuda
julia runbenchmark.jl current cuda

# Compare
julia --project=. compare.jl
```

--------------------------------

### Two-Stage Stochastic Programming Model

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Create a two-stage stochastic program with shared first-stage variables and scenario-specific second-stage variables. Uses `TwoStageExaCore` and `EachScenario`.

```julia
using ExaModels

# Create a two-stage core with 5 scenarios
nscen = 5
c = TwoStageExaCore(nscen; concrete = Val(true))

# First-stage (design) variables - shared across all scenarios
c, d = add_var(c, 3; lvar = 0.0, uvar = 10.0)

# Second-stage (recourse) variables - one copy per scenario
c, v = add_var(c, EachScenario(), 4; lvar = 0.0)

# First-stage constraints
c, design_con = add_con(c, d[i] for i in 1:3; lcon = 1.0)

# Second-stage constraints - replicated per scenario
scenario_data = [(i = i, demand = 10.0 + i) for i in 1:4, s in 1:nscen]
c, recourse_con = add_con(c, EachScenario(),
    v[d.i, s] - d.demand for (d, s) in zip(scenario_data, repeat(1:nscen, inner=4));
    lcon = 0.0
)

# First-stage objective
@add_obj(c, d[i]^2 for i in 1:3)

# Second-stage objective contributions
@add_obj(c, 0.1 * v[i, s]^2 for i in 1:4, s in 1:nscen)

m = ExaModel(c)

# Query scenario structure
println("Number of scenarios: $(get_nscen(m))")
println("Variable scenario indices: $(get_var_scen(m))")
println("Constraint scenario indices: $(get_con_scen(m))")
```

--------------------------------

### ExaModels v0.10 Macro API Usage

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/upgrade.md

Shows the concise macro style for model building with the immutable ExaCore in v0.10. Ensure ExaModels is imported.

```julia
using ExaModels

c = ExaCore(concrete = Val(true))
@add_var(c, x, 10; lvar = 0.0)
@add_obj(c, x[i]^2 for i in 1:10)
m = ExaModel(c)
```

--------------------------------

### Restructuring Dynamics as Data with add_con!

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/constraint_augmentation.md

This method restructures dynamics as data, using `add_con!` to create a single base constraint and an augmentation pattern for all terms, significantly reducing expression patterns and improving GPU performance.

```julia
# Build a single array encoding ALL terms across ALL state equations
terms = [
    (con = state_idx, coeff = coefficient, xi = x_index,
     ui = u_index1, uj = u_index2, vi = v_index)
    for (state_idx, coefficient, x_index, ...) in all_dynamics_terms
]

# One base constraint for all state equations
@add_con(c, dynamics, x[i,j+1,k] - x[i,j,k] for (i,j,k) in state_grid)

# One augmentation pattern for ALL terms across ALL equations
@add_con!(c, dynamics[t.con] += t.coeff * x[t.xi] * u[t.ui] * u[t.uj] * cos(v[t.vi])
    for t in terms
)
```

--------------------------------

### Correct Function: Returning ExaModel

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/patterns.md

Shows the correct way to build a model within a function by returning the final ExaModel, which is constructed from the complete core.

```julia
# CORRECT --- return the ExaModel (or core) to the caller
function build_model()
    c = ExaCore(concrete = Val(true))
    @add_var(c, x, 10)
    @add_obj(c, x[i]^2 for i in 1:10)
    @add_con(c, x[i] + x[i+1] for i in 1:9)
    return ExaModel(c)   # <-- pass the final core to ExaModel
```

--------------------------------

### Add Constraints with += Form

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Demonstrates the use of the += form for adding constraints, which is equivalent to the standard form but more readable. This is useful for accumulating constraints on specific buses.

```julia
gen_data = [(i = i, bus = mod1(i+3, 9)) for i in 1:3]
@add_con!(c, balance[g.bus] += -pg[g.i] for g in gen_data)
```

--------------------------------

### Project.toml Configuration for ExaModels Extension

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/develop.md

This TOML file defines the metadata for an ExaModels extension package, including its name, UUID, authors, dependencies, and testing targets. Ensure ExaModels is listed in [deps].

```toml
# Project.toml

name = "LuksanVlcekModels"
uuid = "0c5951a0-f777-487f-ad29-fac2b9a21bf1"
authors = ["Sungho Shin <sshin@anl.gov>"]
version = "0.1.0"

[deps]
ExaModels = "1037b233-b668-4ce9-9b63-f9f681f55dd2"

[extras]
NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[targets]
test = ["Test", "NLPModelsIpopt"]
```

--------------------------------

### Add Variables to ExaCore Model

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Define optimization variables with dimensions, bounds, and initial values using `add_var` or `@add_var`.

```julia
using ExaModels

c = ExaCore(concrete = Val(true))

# Add 10 variables with default bounds (-Inf, Inf)
c, x = add_var(c, 10)

# Add variables with bounds and initial guess
c, y = add_var(c, 10;
    start = (sin(i) for i = 1:10),
    lvar = 0.0,
    uvar = 10.0
)

# Multi-dimensional variables (9x3 matrix)
c, z = add_var(c, 2:10, 3:5;
    lvar = zeros(9, 3),
    uvar = ones(9, 3)
)

# Using macro syntax (updates core in-place, binds variable name)
c = ExaCore(concrete = Val(true))
@add_var(c, x, 10; lvar = -1, uvar = 1)
@add_var(c, y, 1:5; start = 1.0)

# Access via core.x after macro registration
println(c.x)  # Variable x ∈ R^{10}
```

--------------------------------

### LuksanVlcekModels.jl Extension Module

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/develop.md

This Julia module defines the objective, constraint, and initial guess functions for the Luksan-Vlcek problem, and constructs an ExaModel. It requires ExaModels and exports the model creation function.

```julia
# src/LuksanVlcekModels.jl

module LuksanVlcekModels

import ExaModels

function luksan_vlcek_obj(x,i)
    return 100*(x[i-1]^2-x[i])^2+(x[i-1]-1)^2
end

function luksan_vlcek_con(x,i)
    return 3x[i+1]^3+2*x[i+2]-5+sin(x[i+1]-x[i+2])+sin(x[i+1]+x[i+2])+4x[i+1]-x[i]exp(x[i]-x[i+1])-3
end

function luksan_vlcek_x0(i)
    return mod(i,2)==1 ? -1.2 : 1.0
end

function luksan_vlcek_model(N; backend = nothing)
    
    c = ExaModels.ExaCore(backend, concrete = Val(true))
    c, x = ExaModels.add_variable(
        c, N;
        start = (luksan_vlcek_x0(i) for i=1:N)
    )
    c, _ = ExaModels.add_constraint(
        c,
        luksan_vlcek_con(x,i)
        for i in 1:N-2)
    c, _ = ExaModels.add_objective(c, luksan_vlcek_obj(x,i) for i in 2:N)
    
    return ExaModels.ExaModel(c) # returns the model
end

export luksan_vlcek_model

end # module LuksanVlcekModels
```

--------------------------------

### Add Parameters from Array

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Demonstrates adding parameters to the model from a Julia array. These parameters can be updated later without rebuilding the entire model.

```julia
using ExaModels

c = ExaCore(concrete = Val(true))

# Add parameters from an array
c, theta = add_par(c, [100.0, 1.0, 0.5])
```

--------------------------------

### Macro Expansion for Rebinding

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/patterns.md

Shows how `@add_var`-family macros expand to rebind the core variable in the calling scope.

```julia
c = ExaCore(concrete = Val(true))
@add_var(c, x, 10)          # expands to: c, x = add_var(c, 10)
@add_obj(c, x[i]^2 for i in 1:10)  # expands to: c, _ = add_obj(c, ...)
```

--------------------------------

### GPU Acceleration with KernelAbstractions

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Define a model for NVIDIA GPU execution using KernelAbstractions.jl. Requires CUDA.jl and MadNLPGPU.jl. Note: This code is commented out.

```julia
using ExaModels, KernelAbstractions

# NVIDIA GPU (requires CUDA.jl and MadNLPGPU.jl)
# using CUDA, MadNLPGPU
#
# function create_model_gpu(N)
#     c = ExaCore(Float64; backend = CUDABackend(), concrete = Val(true))
#     d1 = CuArray(1:(N-2))
#     d2 = CuArray(2:N)
#     d3 = CuArray([mod(i,2) == 1 ? -1.2 : 1.0 for i = 1:N])
#
#     @add_var(c, x, N; start = d3)
#     @add_obj(c, 100*(x[i-1]^2 - x[i])^2 + (x[i-1] - 1)^2 for i in d2)
#     @add_con(c, 3x[i+1]^3 + 2*x[i+2] - 5 + 4x[i+1] - x[i]*exp(x[i]-x[i+1]) - 3
#              for i in d1)
#     return ExaModel(c)
# end
#
# m_gpu = create_model_gpu(10000)
# result = madnlp(m_gpu)
```

--------------------------------

### Extract Bound Multipliers

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Explains how to retrieve the dual values associated with the lower and upper bounds of variables using `multipliers_L` and `multipliers_U`. It also shows how to identify active bounds.

```julia
# Get bound multipliers
x_lower = multipliers_L(result, x)  # Lower bound duals
x_upper = multipliers_U(result, x)  # Upper bound duals

# Check which bounds are active
active_lower = findall(x -> x > 1e-6, x_lower)
active_upper = findall(x -> x > 1e-6, x_upper)
```

--------------------------------

### Update Parameters Without Rebuilding Model

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Demonstrates how to update the values of existing parameters after the model has been created but before it's solved. This is crucial for scenario analysis or warm-starting.

```julia
# Update parameters without rebuilding model
m = ExaModel(c)
set_parameter!(c, penalty, [300.0, 75.0])
```

--------------------------------

### Solve Luksan-Vlcek Problem with ExaModels.jl

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Defines and solves the Luksan-Vlcek nonlinear programming test problem with a specified number of variables. Requires ExaModels and NLPModelsIpopt packages.

```julia
using ExaModels, NLPModelsIpopt

function solve_luksan_vlcek(N)
    c = ExaCore(concrete = Val(true))

    # Variables with alternating initial guess
    @add_var(c, x, N; start = (mod(i, 2) == 1 ? -1.2 : 1.0 for i = 1:N))

    # Rosenbrock-like objective
    @add_obj(c, 100 * (x[i-1]^2 - x[i])^2 + (x[i-1] - 1)^2 for i = 2:N)

    # Nonlinear equality constraints
    @add_con(c,
        3x[i+1]^3 + 2*x[i+2] - 5 +
        sin(x[i+1] - x[i+2]) * sin(x[i+1] + x[i+2]) +
        4x[i+1] - x[i] * exp(x[i] - x[i+1]) - 3
        for i = 1:(N-2)
    )

    m = ExaModel(c)
    result = ipopt(m; print_level = 5, max_iter = 100)

    return (
        status = result.status,
        objective = result.objective,
        iterations = result.iter,
        solution = solution(result, x)
    )
end

# Solve with N=100 variables
result = solve_luksan_vlcek(100)
println("Status: $(result.status)")
println("Final objective: $(result.objective)")
println("Iterations: $(result.iterations)")
```

--------------------------------

### Add a Parameter to ExaCore

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/parameters.md

Add a parameter to the ExaCore by providing its name and a vector of initial values. Parameters can be used in objective functions and constraints.

```julia
@add_parameter(c_param, θ, [100.0, 1.0])  # [penalty_coeff, offset]
```

--------------------------------

### Register and Use Custom Bivariate Function

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Register a custom bivariate function and its derivatives for use in an ExaModel. Ensure all required derivatives are provided.

```julia
softplus_sum(x, y) = log(1 + exp(x + y))
dsoftplus_sum_x(x, y) = exp(x + y) / (1 + exp(x + y))
dsoftplus_sum_y(x, y) = exp(x + y) / (1 + exp(x + y))
ddsoftplus_sum_xx(x, y) = exp(x + y) / (1 + exp(x + y))^2
ddsoftplus_sum_xy(x, y) = exp(x + y) / (1 + exp(x + y))^2
ddsoftplus_sum_yy(x, y) = exp(x + y) / (1 + exp(x + y))^2
@register_bivariate(softplus_sum, dsoftplus_sum_x, dsoftplus_sum_y,
                    ddsoftplus_sum_xx, ddsoftplus_sum_xy, ddsoftplus_sum_yy)

c = ExaCore(concrete = Val(true))
@add_var(c, x, 5; start = 1.0)
@add_var(c, y, 5; start = 1.0)
@add_obj(c, relu3(x[i]) + softplus_sum(x[i], y[i]) for i = 1:5)
m = ExaModel(c)
```

--------------------------------

### Modify Penalty Coefficient and Re-solve

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/parameters.md

This code snippet demonstrates how to change the penalty coefficient for the IPOPT solver and then re-run the optimization. It prints the objective value of the new solution.

```julia
set_parameter!(c_param, θ, [200.0, 0.5])  # Change the offset in the objective
result3 = ipopt(m_param)
println("Modified offset objective: $(result3.objective)")
```

--------------------------------

### Extract Primal Solution for Variables

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Shows how to retrieve the optimal values for the model's variables after a solution has been found. This uses the `solution` function with the variable object.

```julia
using ExaModels, NLPModelsIpopt

c = ExaCore(concrete = Val(true))
@add_var(c, x, 1:10; lvar = -1, uvar = 1)
@add_obj(c, (x[i] - 2)^2 for i = 1:10)
@add_con(c, g, x[i] + x[i+1] for i = 1:9; lcon = -1, ucon = 1)

m = ExaModel(c)
result = ipopt(m; print_level = 0)

# Get primal solution for variables
x_sol = solution(result, x)
println("Solution x: $x_sol")
```

--------------------------------

### Register Custom Univariate Function

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Shows how to register a custom univariate function (relu3) along with its first and second derivatives using `@register_univariate`. This allows the function to be used within ExaModels expressions.

```julia
using ExaModels

# Register a custom univariate function
relu3(x) = x > 0 ? x^3 : zero(x)
drelu3(x) = x > 0 ? 3*x^2 : zero(x)
ddrelu3(x) = x > 0 ? 6*x : zero(x)
@register_univariate(relu3, drelu3, ddrelu3)
```

--------------------------------

### Add Constraints to ExaCore Model

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Define constraints using `add_con` or `@add_con` with specified bounds. Supports equality and inequality constraints.

```julia
using ExaModels

c = ExaCore(concrete = Val(true))
@add_var(c, x, 10)

# Equality constraints (default lcon=ucon=0)
c, con = add_con(c, x[i] + x[i+1] for i = 1:9)

# Inequality constraints with bounds
c, con2 = add_con(c, x[i]^2 for i = 1:10;
    lcon = -1.0,
    ucon = (i * 0.5 for i = 1:10)
)

# Upper-bounded only (lower = -Inf)
@add_con(c, x[i]^2 for i = 1:10; lcon = -Inf, ucon = 1.0)

# Named constraints accessible via core.name
@add_con(c, g, x[i] - x[i+1] for i = 1:9; lcon = -0.5, ucon = 0.5)
println(c.g)  # Constraint g

# Create empty constraint slots for later augmentation
c = ExaCore(concrete = Val(true))
@add_var(c, x, 10)
c, balance = add_con(c, 9; lcon = -1.0, ucon = 1.0)
```

--------------------------------

### add_con! Syntax: Two-Argument += Form

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/constraint_augmentation.md

The two-argument `+=` form of `add_con!` embeds the constraint and index directly, offering a more natural reading. This syntax is fully equivalent to the three-argument form.

```julia
@add_con!(c, g[a.bus] += p[a.i] for a in arc_data)
```

```julia
c, _ = add_con!(c, g[a.bus] += p[a.i] for a in arc_data)
```

--------------------------------

### Traditional Dynamic Expression Pattern

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/constraint_augmentation.md

This approach creates a distinct expression pattern for each dynamic term, leading to slow compilation and poor GPU performance due to numerous kernel launches.

```julia
for i in 1:20
    @add_con(c, x[i,j+1,k] - x[i,j,k] - dyns[i](x,u,v,j+1,k) * dt
             for (j,k) in grid)
end
```

--------------------------------

### Use Parameters in Objective Function

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Defines variables and then uses previously added parameters within the objective function. This allows for flexible model definition where coefficients can be changed.

```julia
# Use parameters in objective
@add_var(c, x, 10; start = 1.0)
@add_obj(c, penalty[1] * (x[i]^2 - 1)^2 + penalty[2] * x[i] for i = 1:10)
```

--------------------------------

### Incorrect Function: Forgetting to Return Core

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/patterns.md

Illustrates a common pitfall where a function modifies a local core but does not return it, leading to data loss for the caller.

```julia
# WRONG --- the caller gets nothing useful
function build_model_wrong()
    c = ExaCore(concrete = Val(true))
    @add_var(c, x, 10)
    @add_obj(c, x[i]^2 for i in 1:10)
    @add_con(c, x[i] + x[i+1] for i in 1:9)
    # c is local --- it is lost when the function returns!
```

--------------------------------

### Augment Constraints in ExaCore Model

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Add expression terms to existing constraint rows using `add_con!` or `@add_con!`. Useful for GPU performance.

```julia
using ExaModels

c = ExaCore(concrete = Val(true))
@add_var(c, x, 10)
@add_var(c, p, 5)
@add_var(c, pg, 3)

# Create base constraint
@add_con(c, balance, x[i] for i = 1:9; lcon = 0.0, ucon = 0.0)

# Three-argument form: add terms to specific constraint rows
arc_data = [(i = i, bus = mod1(i, 9)) for i in 1:5]
@add_con!(c, balance, a.bus => p[a.i] for a in arc_data)
```

--------------------------------

### Define Objective Function with Parameters

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/parameters.md

Incorporate parameters into the objective function definition. Here, θ[1] and θ[2] are used as coefficients and offsets.

```julia
@add_objective(c_param, θ[1] * (x_p[i-1]^2 - x_p[i])^2 + (x_p[i-1] - θ[2])^2 for i = 2:N)
```

--------------------------------

### Add Variables to ExaCore

Source: https://github.com/exanauts/examodels.jl/blob/main/docs/src/parameters.md

Define variables for the ExaCore, specifying their dimension and initial values. These variables can interact with parameters in the model.

```julia
N = 10
@add_variable(c_param, x_p, N; start = (mod(i, 2) == 1 ? -1.2 : 1.0 for i = 1:N))
```

--------------------------------

### Use Subexpression in Objective and Constraints

Source: https://context7.com/exanauts/examodels.jl/llms.txt

Utilizes the previously defined subexpression 's' in both the objective function and constraints. This demonstrates how subexpressions can simplify complex model formulations.

```julia
# Use subexpression in objective and constraints
@add_obj(c, (s[i] - 1)^2 for i in 1:N)
@add_con(c, s[i] + s[i+1] for i in 1:(N-1); lcon = 0.0)
```