# Minimisation with equality constraints

**URL:** <https://discourse.julialang.org/t/minimisation-with-equality-constraints/107358>\
**Category:** Optimization (Mathematical)\
**Tags:** question, package, syntax, optimization\
**Created:** [December 9, 2023, 1:40pm UTC](https://discourse.julialang.org/t/minimisation-with-equality-constraints/107358 "2023-12-09T13:40:40Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![eveningsilverfox](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eveningsilverfox/32/51932_2.png) [@eveningsilverfox](https://discourse.julialang.org/u/eveningsilverfox)\
**Post date:** [December 9, 2023, 1:40pm UTC](https://discourse.julialang.org/t/minimisation-with-equality-constraints/107358/1 "2023-12-09T13:40:40Z")

</div>

I am trying a constrained minimisation with equality constraints.

```julia
using LinearAlgebra
using Optimization, OptimizationMOI, AmplNLWriter, Ipopt_jll

function objectivehelp(N)
    car = rand(N) .+ 1; sar = rand(N);

    return car, sar
end

function objective(x, N)
    car, sar = objectivehelp(N);
    res = 0;
    for ij = 1:N
        res = res + car[ij] * (x[ij]-sar[ij])^2;
    end
    res = abs(res);

    return res
end

function constraint(x,N)
    # res = 0;
    # for ij = 1:N
    # res = res + x[ij]^2;
    # end
    res = zeros(Float64, 1);
    for ij = 1:N
        res[1] = res[1] + x[ij]^2;
    end
    
    return res
end

N = 3;
car = rand(N) .+ 1; sar = rand(N);

obj(x, p) = objective(x,N)
cons(res, x, p) = (res .= constraint(x,N))

x0 = ones(Float64, N);

optprob = OptimizationFunction(obj, Optimization.AutoModelingToolkit(), cons = cons)
prob = OptimizationProblem(optprob, x0, lcons = ones(Float64,1), ucons = ones(Float64,1))
sol = solve(prob, AmplNLWriter.Optimizer(Ipopt_jll.amplexe), maxiters = 1000, abstol = 1e-6, reltol = 1e-6)

```

The code using the commented part in the constraint function works fine, but when I use the uncommented part instead, it throws the error

```julia
ERROR: LoadError: MethodError: no method matching Float64(::Symbolics.Num)

Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
   @ Base rounding.jl:207
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:792
  (::Type{T})(::SymbolicUtils.Symbolic) where T<:Union{AbstractFloat, Integer, Complex{<:AbstractFloat}, Complex{<:Integer}}
   @ Symbolics ~/.julia/packages/Symbolics/5BibL/src/Symbolics.jl:149
  ...

Stacktrace:
  [1] convert(#unused#::Type{Float64}, x::Symbolics.Num)
    @ Base ./number.jl:7
  [2] setindex!(A::Vector{Float64}, x::Symbolics.Num, i1::Int64)
    @ Base ./array.jl:969
  [3] constraint(x::Vector{Symbolics.Num}, N::Int64)
    @ Main ~/Work/Higgs work and codes/Codes/testoptim.jl:28
  [4] cons(res::Vector{Symbolics.Num}, x::Vector{Symbolics.Num}, p::Vector{Any})
    @ Main ~/Work/Higgs work and codes/Codes/testoptim.jl:38
  [5] modelingtoolkitize(prob::OptimizationProblem{true, OptimizationFunction{true, AutoModelingToolkit, typeof(obj), Nothing, Nothing, Nothing, typeof(cons), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float64}, SciMLBase.NullParameters, Nothing, Nothing, Nothing, Vector{Float64}, Vector{Float64}, Nothing, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/dCa81/src/systems/optimization/modelingtoolkitize.jl:23
  [6] modelingtoolkitize(prob::OptimizationProblem{true, OptimizationFunction{true, AutoModelingToolkit, typeof(obj), Nothing, Nothing, Nothing, typeof(cons), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float64}, SciMLBase.NullParameters, Nothing, Nothing, Nothing, Vector{Float64}, Vector{Float64}, Nothing, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/dCa81/src/systems/optimization/modelingtoolkitize.jl:6
  [7] instantiate_function(f::Function, cache::Optimization.ReInitCache{Vector{Float64}, SciMLBase.NullParameters}, adtype::AutoModelingToolkit, num_cons::Int64)
    @ OptimizationMTKExt ~/.julia/packages/Optimization/fPVIW/ext/OptimizationMTKExt.jl:58
  [8] OptimizationMOI.MOIOptimizationNLPCache(prob::OptimizationProblem{true, OptimizationFunction{true, AutoModelingToolkit, typeof(obj), Nothing, Nothing, Nothing, typeof(cons), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float64}, SciMLBase.NullParameters, Nothing, Nothing, Nothing, Vector{Float64}, Vector{Float64}, Nothing, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, opt::AmplNLWriter.Optimizer; callback::Nothing, kwargs::Base.Pairs{Symbol, Union{Nothing, Real}, NTuple{4, Symbol}, NamedTuple{(:maxiters, :maxtime, :abstol, :reltol), Tuple{Int64, Nothing, Float64, Float64}}})
    @ OptimizationMOI ~/.julia/packages/OptimizationMOI/Ip1EA/src/nlp.jl:109
  [9] MOIOptimizationNLPCache
    @ ~/.julia/packages/OptimizationMOI/Ip1EA/src/nlp.jl:105 [inlined]
 [10] __init(prob::OptimizationProblem{true, OptimizationFunction{true, AutoModelingToolkit, typeof(obj), Nothing, Nothing, Nothing, typeof(cons), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float64}, SciMLBase.NullParameters, Nothing, Nothing, Nothing, Vector{Float64}, Vector{Float64}, Nothing, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, opt::AmplNLWriter.Optimizer; maxiters::Int64, maxtime::Nothing, abstol::Float64, reltol::Float64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ OptimizationMOI ~/.julia/packages/OptimizationMOI/Ip1EA/src/OptimizationMOI.jl:180
 [11] __init
    @ ~/.julia/packages/OptimizationMOI/Ip1EA/src/OptimizationMOI.jl:172 [inlined]
 [12] #init#617
    @ ~/.julia/packages/SciMLBase/ynHlA/src/solve.jl:163 [inlined]
 [13] init
    @ ~/.julia/packages/SciMLBase/ynHlA/src/solve.jl:161 [inlined]
 [14] solve(::OptimizationProblem{true, OptimizationFunction{true, AutoModelingToolkit, typeof(obj), Nothing, Nothing, Nothing, typeof(cons), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float64}, SciMLBase.NullParameters, Nothing, Nothing, Nothing, Vector{Float64}, Vector{Float64}, Nothing, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::AmplNLWriter.Optimizer; kwargs::Base.Pairs{Symbol, Real, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:maxiters, :abstol, :reltol), Tuple{Int64, Float64, Float64}}})
    @ SciMLBase ~/.julia/packages/SciMLBase/ynHlA/src/solve.jl:94
 [15] top-level scope
    @ ~/Work/Higgs work and codes/Codes/testoptim.jl:44
 [16] include(fname::String)
    @ Base.MainInclude ./client.jl:478
 [17] top-level scope
    @ REPL[1]:1
in expression starting at /home/aritra/Work/Higgs work and codes/Codes/testoptim.jl:44

```

I was looking for a way to supply multiple equality constraints, and hence I was trying out the commented part.

Any ideas what the error is?

---

<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:** [December 9, 2023, 8:37pm UTC](https://discourse.julialang.org/t/minimisation-with-equality-constraints/107358/2 "2023-12-09T20:37:00Z")

</div>

Have you considered using JuMP for this?

```julia
using JuMP, Ipopt
N = 3
car = rand(N) .+ 1
sar = rand(N)
model = Model(Ipopt.Optimizer)
@variable(model, x[1:N])
@objective(model, Min, sum(car[i] * (x[i] - sar[i])^2 for i in 1:N))
@constraint(model, sum(x[i]^2 for i in 1:N) == 1)

```

---

<div class="post-metadata">

**Author:** ![eveningsilverfox](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eveningsilverfox/32/51932_2.png) [@eveningsilverfox](https://discourse.julialang.org/u/eveningsilverfox)\
**Post date:** [December 9, 2023, 10:45pm UTC](https://discourse.julialang.org/t/minimisation-with-equality-constraints/107358/3 "2023-12-09T22:45:25Z")

</div>

> [@odow](#):
>
> ```julia
> N = 3
> car = rand(N) .+ 1
> sar = rand(N)
> model = Model(Ipopt.Optimizer)
> @variable(model, x[1:N])
> @objective(model, Min, sum(car[i] * (x[i] - sar[i])^2 for i in 1:N))
> @constraint(model, sum(x[i]^2 for i in 1:N) == 1)
> 
> ```

Thanks. Is it possible to minimise user defined functions as shown below? The example provided below doesn’t work. I have to minimise such explicitly-stated functions, which cannot be written in one line.

```julia
using LinearAlgebra
using JuMP, Ipopt

function objectivehelp(x, N)
    xp = zeros(N);
    for ij = 1:N
        xp[ij] = 2*x[ij];
    end

    return xp
end

function objective(x, N)
    car = rand(N) .+ 1; sar = rand(N);
    xp = objectivehelp(x, N);
    
    res = 0;
    for ij = 1:N
        res = res + car[ij] * ( xp[ij]^0.5-sar[ij])^2;
    end
    res = abs(res);

    return res
end

function constrainthelp(x)
    res = 2*x;

    return res
end

function constraint(x,N)
    y = constrainthelp(x);
    cmat = zeros(Float64, 2,N);
    for ij = 1:2
        for jk = 1:N
            cmat[ij,jk] = rand()-0.5;
        end
    end
    res = cmat * y;
    
    return res
end

obj(x) = objective(x,N)
cons(x) = constraint(x,N)

N = 3
car = rand(N) .+ 1
sar = rand(N)
model = Model(Ipopt.Optimizer)
@variable(model, x[1:N])
@objective(model, Min, obj(x))
@constraint(model, cons(x))
optimize!(model)
objective_value(model)
value(x)

```

---

<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:** [December 10, 2023, 1:07am UTC](https://discourse.julialang.org/t/minimisation-with-equality-constraints/107358/4 "2023-12-10T01:07:50Z")

</div>

A few comments:

- Your functions need to be deterministic. Currently, each time you call `objective` and `constraint` you will get a different answer
- `zeros(N)` creates a `Vector{Float64}`. You can use Julia’s type dispatch to create a `zero(T, N)` instead, which creates an array of the right type.
- Your `@constraint` line is missing a comparison operator like `== 1`.
- Your functions include `x^0.5`, for which the gradient isn’t defined when `x = 0`. You probably need to set a non-zero lower bound for `x`

Making all those changes, I get something like this:

```julia
julia> using JuMP, Ipopt

julia> function objective_help(x::Vector{VariableRef})
           N = length(x)
           xp = zeros(AffExpr, N)
           for ij in 1:N
               xp[ij] = 2 * x[ij]
           end
           return xp
       end
objective_help (generic function with 2 methods)

julia> function objective(x, car, sar)
           N = length(x)
           xp = objective_help(x)
           res = sum(car[ij] * (xp[ij]^0.5 - sar[ij])^2 for ij in 1:N)
           return abs(res)
       end
objective (generic function with 1 method)

julia> constraint_help(x) = 2 .* x
constraint_help (generic function with 1 method)

julia> function constraint(x, cmat)
           y = constraint_help(x);
           return cmat * y
       end
constraint (generic function with 1 method)

julia> N = 3
3

julia> car = rand(N) .+ 1
3-element Vector{Float64}:
 1.3146682360894872
 1.0688322590056367
 1.213419895106776

julia> sar = rand(N)
3-element Vector{Float64}:
 0.22311572604394314
 0.525736516469558
 0.0945780732842676

julia> cmat = zeros(Float64, 2, N)
2×3 Matrix{Float64}:
 0.0 0.0 0.0
 0.0 0.0 0.0

julia> for ij in 1:2, jk in 1:N
           cmat[ij, jk] = rand() - 0.5
       end

julia> model = Model(Ipopt.Optimizer)
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Ipopt

julia> @variable(model, x[1:N] >= 0.00001)
3-element Vector{VariableRef}:
 x[1]
 x[2]
 x[3]

julia> @objective(model, Min, objective(x, car, sar))
abs((1.3146682360894872 * ((((2 x[1]) ^ 0.5) - 0.22311572604394314) ^ 2.0)) + (1.213419895106776 * ((((2 x[3]) ^ 0.5) - 0.0945780732842676) ^ 2.0)) + (1.0688322590056367 * ((((2 x[2]) ^ 0.5) - 0.525736516469558) ^ 2.0)))

julia> @constraint(model, constraint(x, cmat) .== 1)
2-element Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape}}:
 -0.018778023797837218 x[1] + 0.7589910154986776 x[2] - 0.8148125363688683 x[3] = 1
 -0.6298413441703727 x[1] - 0.17177726380301261 x[2] + 0.9023777745930739 x[3] = 1

julia> optimize!(model)
This is Ipopt version 3.14.13, running with linear solver MUMPS 5.6.1.

Number of nonzeros in equality constraint Jacobian...: 6
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 6

Total number of variables............................: 3
                     variables with only lower bounds: 3
                variables with lower and upper bounds: 0
                     variables with only upper bounds: 0
Total number of equality constraints.................: 2
Total number of inequality constraints...............: 0
        inequality constraints with only lower bounds: 0
   inequality constraints with lower and upper bounds: 0
        inequality constraints with only upper bounds: 0

iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
   0 1.6923605e-01 1.00e+00 2.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
   1 1.8237829e-01 9.95e-01 1.58e+01 -1.0 1.81e+00 - 7.62e-03 5.47e-03h 1
   2 7.9818312e+00 2.22e-16 3.72e+02 -1.0 3.14e+00 - 6.50e-03 1.00e+00h 1
   3 8.0084526e+00 4.44e-16 8.85e-01 -1.0 1.09e-02 - 1.00e+00 1.00e+00f 1
   4 8.0059375e+00 4.44e-16 1.01e-03 -1.7 8.00e-04 - 1.00e+00 1.00e+00f 1
   5 7.9808356e+00 2.22e-16 4.95e-01 -3.8 1.11e-02 - 1.00e+00 1.00e+00f 1
   6 7.9803142e+00 0.00e+00 2.33e-02 -3.8 8.93e-04 - 1.00e+00 1.00e+00f 1
   7 7.9802377e+00 2.22e-16 5.30e-03 -3.8 3.58e-04 - 1.00e+00 1.00e+00f 1
   8 7.9802365e+00 1.11e-16 7.93e-06 -3.8 1.30e-05 - 1.00e+00 1.00e+00f 1
   9 7.9802314e+00 2.22e-16 6.59e-04 -5.7 1.17e-04 - 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
  10 7.9802314e+00 2.22e-16 2.93e-08 -5.7 7.60e-07 - 1.00e+00 1.00e+00h 1
  11 7.9802314e+00 2.22e-16 1.09e-07 -8.6 1.46e-06 - 1.00e+00 1.00e+00h 1
  12 7.9802314e+00 1.11e-16 2.60e-14 -8.6 1.23e-10 - 1.00e+00 1.00e+00h 1

Number of Iterations....: 12

                                   (scaled) (unscaled)
Objective...............: 3.2028226115109715e+00 7.9802314272782713e+00
Dual infeasibility......: 2.6014114945468764e-14 6.4817407275241498e-14
Constraint violation....: 1.1102230246251565e-16 1.1102230246251565e-16
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 2.5059035775423343e-09 6.2437708574182410e-09
Overall NLP error.......: 2.5059035775423343e-09 6.2437708574182410e-09

Number of objective function evaluations = 13
Number of objective gradient evaluations = 13
Number of equality constraint evaluations = 13
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 1
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 12
Total seconds in IPOPT = 0.006

EXIT: Optimal Solution Found.

julia> objective_value(model)
7.980231427278271

julia> value.(x)
3-element Vector{Float64}:
 0.004325149667465344
 3.1554184142833086
 1.7118698431331894

```

In many cases you really should consider writing out the functions algebraically though. if only to think carefully about what you are actually modelling.

---

<div class="post-metadata">

**Author:** ![eveningsilverfox](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eveningsilverfox/32/51932_2.png) [@eveningsilverfox](https://discourse.julialang.org/u/eveningsilverfox)\
**Post date:** [December 10, 2023, 2:47pm UTC](https://discourse.julialang.org/t/minimisation-with-equality-constraints/107358/5 "2023-12-10T14:47:19Z")

</div>

> [@odow](#):
>
> `xp = zeros(AffExpr, N)`

Thanks.  
My actual functions are deterministic, the random coefficients car and sar were created just for easy instantiation.

So the data types in the “intermediate” functions have to be AffExpr.  
Is there a way of getting around this? At the moment I have a very complicated function (defined in terms of other functions), and it seems that to make it work with JuMP I have to restructure the entire code from scratch.  
From your answer, it seems that I can use the data type T instead, right?

---

<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:** [December 10, 2023, 9:13pm UTC](https://discourse.julialang.org/t/minimisation-with-equality-constraints/107358/6 "2023-12-10T21:13:06Z")

</div>

> So the data types in the “intermediate” functions have to be AffExpr.

No. This is just an example. For `zeros(T, N)`, you should use the appropriate `T`, whatever that may be.

You could also do:

```julia
xp = [2 * x[ij] for ij in 1:N]

```

and let Julia figure out the correct type. You’re only having a problem because you are trying to pre-allocate the buffer `xp`.

The pre-allocation is also the cause of your original issue. You used  
`res = zeros(Float64, 1);` but `x[ij]^2` is not a `Float64`, it was a `Symbolics.Num`.

> At the moment I have a very complicated function (defined in terms of other functions), and it seems that to make it work with JuMP I have to restructure the entire code from scratch.

You could also register a user-defined operator:

> **[Nonlinear Modeling · JuMP](https://jump.dev/JuMP.jl/stable/manual/nonlinear/#jump_user_defined_operators)**
>
> Documentation for JuMP.

But if your function is very complicated, see:

> **[Should you use JuMP? · JuMP](https://jump.dev/JuMP.jl/stable/should_i_use/#You-want-to-optimize-a-complicated-Julia-function)**
>
> Documentation for JuMP.

What you should use depends on whether you can automatically differentiate the function, so if you want advice, perhaps you could post a more realistic example? It’s hard to offer general advised based on the simple example above.